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Abstract 

The next-to-next- to-leading order post-Newtonian spin-orbit and spin(l)-spin(2) Hamiltonians for bi- 
nary compact objects in general relativity are derived. The Arnowitt-Deser-Misner canonical formalism 
and its generalization to spinning compact objects in general relativity are presented and a fully reduced 
matter-only Hamiltonian is obtained. Several simplifications using integrations by parts are discussed. 
Approximate solutions to the constraints and evolution equations of motion are provided. Technical 
details of the integration procedures are given including an analysis of the short-range behavior of the 
integrands around the sources. The Hamiltonian of a test-spin moving in a stationary Kerr spacetime is 
obtained by rather simple approach and used to check parts of the mentioned results. Kinematical consis- 
tency checks by using the global (post-Newtonian approximate) Poincare algebra are applied. Along the 
way a self-contained overview for the computation of the 3PN ADM point-mass Hamiltonian is provided, 
too. 

PACS numbers: 04.25.Nx, 04.20.Fy, 04.25.-g, 97.80.-d, 45.50.Jf 

Keywords: post-Newtonian approximation; canonical formalism; approximation methods; equations of 
motion; binary stars 



1 Introduction 

Since Einstein discovered the theory of general relativity [1-5], many attempts to solve the field equations 
were undertaken. Yet only a few analytical solutions to the full field equations are known at the time being, 
mostly for highly symmetric matter or field configurations [6, 7]. The most famous solutions of the field 
equations are the Schwarzschild [8] and the Kerr [9] ones. Even so, an analytic solution for binary systems 
of black holes (or even neutron stars) is missing and unlikely to be found in the future. Nevertheless, such 
binaries are very interesting. In particular they constitute the most promising and strongest sources for 
gravitational waves, one of the most fascinating predictions of the theory of general relativity [10, 11]. 

To observe gravitational waves one needs very sensitive detectors because of the tiny cross section of 
the waves with matter There exist several ground-based detector projects like e.g. Geo600, VIRGO, and 
LIGO [12-14] for this purpose. Their sensitivity increased during the last years due to continuous upgrades 
and they probably will detect gravitational waves directly within the next/ew years. The mentioned large- 
scale detectors were started to be built after gravitational waves had been observed indirectly by measuring 
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the change of the radial orbital period e.g. for the binary pulsar PSR 1913+16 by Hulse and Taylor [15] in 
1978 (Nobel prize 1993). Further precise observations of orbital period decays are the double pulsar PSR 
J0737-3039, [16, 17], and the white dwarf binary J0651+2844 [18]. Their theoretical prediction is based 
on the famous quadrupole formula (see e.g. [19, Eq. (4.12)]) which gives an expression for the orbital 
averaged energy loss of the whole system due to gravitational radiation and can thus be translated into a 
decay of the orbital period. 

If gravitational wave signals are detected in the near future, a great effort in data analysis will be 
necessary to extract them from the noisy raw data. Necessary ingredients to achieve this goal on the way 
towards gravitational wave astronomy are predictions on which kind of signals can be expected, e.g. from 
binary sources. The theoretical form of the signals - which should incorporate all known physical effects 
to some specific order of magnitude - is called template and relies on a solution of the field equations 
and therefore on the physical parameters of the system. Unfortunately, after around one hundred years of 
research there are no dynamical analytical solutions of the field equations for an n-body system (n > 2) 
known. There are two possibilities to circumvent this problem. The first one is to rely on numerical 
simulations. The second possibility is to rely on approximation methods to extract solutions for these kind 
of problems from the field equations. These include in particular the post-Minkowskian approximation, the 
post-Newtonian (FN) approximation, extreme mass ratios (including testmass case), and the effective-one- 
body approach. 

One of the most successful approximation methods is the post-Newtonian approximation, a slow motion 
and wide separation approximation. It is used to treat the finite propagation velocity of the gravitational 
field approximately as an instantaneous effect and therefore "freeze" its dynamics. Thus the field degrees 
of freedom are eliminated in this approximation scheme. Afterwards one is left with ordinary differential 
equations for positions, momenta, and spins of the objects in the system only. It is convenient to encode 
these equations of motion into a Lagrangian or a Hamiltonian. 

To get a more quantitative understanding of the post-Newtonian approximation, consider a gravitation- 
ally bound binary system in the Newtonian limit. In this case one can relate the Newtonian kinetic energy 
and the Newtonian gravitational potential through the Virial theorem, namely 

GM Ts 
c c r r 

with V denoting a typical orbital velocity in the system, c meaning the speed of light, G Newton's grav- 
itational constant, M the total mass of the system, r a typical distance between the gravitating objects, 
and Vs being the Schwarzschild radius of the total system. Every order in w^/c^ is denoted as one relative 
post-Newtonian order The post-Newtonian approximation was used to obtain matter equations of motion 
for some special systems akeady shortly after general relativity had been developed, see e.g. the IPN bi- 
nary Lagrangian in [20] and the IPN equations of motion in e.g. [21]. Furthermore previous results for 
non-spinning objects obtained within the formalism used in the present article, namely the canonical for- 
malism of Arnowitt, Deser, and Misner (ADM) [22], are the 2PN [23-25], 2.5PN [19, 26], 3PN [27-31], 
and 3.5PN [32, 33] Hamiltonians. For various other (non-canonical) derivations of post-Newtonian results 
for non-rotating objects see [34^2] and references therein. 

Regarding the spin corrections to the post-Newtonian approximation, the leading order can be found 
in [43-46]. Interestingly the leading-order equations of motion were already obtained earlier within the 
(more general) post-Minkowskian approximation [47, 48]. The post-Newtonian next-to-leading order in 
spin was only tackled more recently [49-58]. The post-Newtonian next-to-next-to-leading order spin-orbit 
and spin(l)-spin(2) Hamiltonians are the subject of the present paper For very rapidly rotating objects they 
can be comparable in strength to 3.5PN and 4PN corrections, respectively. Half a post-Newtonian order 
above them the first (leading-order) spin-dependent radiative Hamiltonians appear The spin-orbit and 
spin(l)-spin(2) ones were already obtained recently [59]. At the 3.5PN level one should further include all 
Hamiltonians cubic in the spins derived in [60, 61]. Notice that these cubic Hamiltonians are only known 
for binary black holes so far, whereas all other mentioned Hamiltonians (including the ones derived in the 
present paper) are valid for general compact objects [or have been generalized to this case, see [62-65] 
for the spin(l)-spin(l) level]. However, the Hamiltonians in the present paper are only obtained for binary 
systems (but many results for three and more objects already exist [66-70]). Besides spin effects, tidal 
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contributions to the post-Newtonian approximation become very important for general compact objects 
like neutron stars [71-73], also see e.g. [74] for the extreme mass ratio case. 

The present article provides a detailed exposition of how the next-to-next-to-leading order (NNLO) 
spin-orbit and spin(l)-spin(2) post-Newtonian Hamiltonians can be derived from an extension of the canon- 
ical formalism of Arnowitt, Deser, and Misner [22] . The mentioned extension of the ADM formalism refers 
to a generalization from (non-rotating) point-masses (PM) to rotating objects [75], see also [76-78]. Re- 
sults for these Hamiltonians were already presented in this journal [79, 80]. Their technicahties are also 
discussed in the (german) PhD thesis of JH [81]. A corresponding Lagrangian potential for the NNLO 
spin(l)-spin(2) interaction was derived simultaneously by M. Levi in [82] via an effective field theory 
(EFT) approach. A comparisons between EFT and ADM results at NNLO spin(l)-spin(2) was not yet 
undertaken, as it is not straightforward (see [83, 84] for a discussion at NLO). However, in [85, 86] the 
equations of motion in harmonic gauge were calculated at NNLO spin-orbit level and agreement with our 
Hamiltonian was found. 

For unbound systems one can in general not relate the velocities of the objects to the strength of the 
gravitational coupling, i.e. the basic relation (1) of the post-Newtonian approximation is not applicable. 
For these kinds of systems (e.g. scattering of black holes) an useful approximation is the so called post- 
Minkowskian approximation, which is an expansion in powers of the gravitational coupling constant G 
only and thus also appropriate for very high velocities and weak fields. The first post-Minkowskian ap- 
proximation for non-rotating objects was used to derive the Hamiltonian in [87] within the ADM formal- 
ism. In principle the expressions given in [26] can be used to derive the ADM Hamiltonian in the second 
post-Minkowskian approximation if there is no incoming radiation. But the integrals have not been given 
in closed form yet and since we are only interested in gravitationally bound systems, we retreat to the 
post-Newtonian approximation. 

The Hamiltonians derived in this article are not the end of the journey. In order to extract useful 
information (i.e. the parameters of a binary) from the gravitational waves, one needs the solution of the 
post-Newtonian equations of motion for the binary. ' For known orbital parameterizations one can further 
calculate the far-zone radiation field (see e.g [98-100] for general formalism of treating radiation in general 
relativity, [101] for higher order radiation losses in point-mass binaries, [102-104] for spin effects on the 
radiation, [105] for spin-dependent tail effects, and [106, 107] for multipole moments including spin up 
to 2.5PN). In case of eccentric orbits the radiation consists of several modes which may be extracted by 
a mode decomposition (see e.g. [108-110] for decomposition of the radiation field in tensor spherical 
harmonics, and computing and solving ordinary differential equations for mean motion and eccentricity in 
a binary without spin; see also [11 1-1 17] for higher order mode decomposition of multipole moments also 
for point-mass binaries only). Further for known parameterizations one is also able to calculate the loss 
of energy and angular momentum during the inspiral process, see e.g. [1 14, 118, 1 19] for time evolution 
effects on the template banks and [59, 120] for the near-zone luminosity. This is necessary to construct the 
mentioned template banks to extract the physical parameters from a noisy signal via a matched filtering 
procedure. There the "fitting factor" is a very sensitive indicator for the performance of the template bank 
vis-a-vis the real signal. An introduction to matched filtering can be found in [121, 122]. There exists a 
plethora of articles referring to circular inspiral without spin, e.g. [123-126]. 

Though the effects considered in the present article are very small, they still probably have a serious im- 
pact on future template banks. The reason for this is that even tiny contributions to the binary's interaction 
accumulate during the long inspiral phase (where the post-Newtonian approximation is still valid) and thus 
may become observable for potentially planned space based detectors in the future. During the very late 
inspiral phase these effects will become more important, but the post-Newtonian approximation will break 
down due to the highly nonlinear behavior of the dynamics and high velocities (f /c > 1 /3). To overcome 
this problem it is most convenient to extrapolate to this nonlinear regime by resumming the post-Newtonian 



'The following literature and in particular [88] gives a complete overview over the research area of parameterization. See e.g. 
[89] for a point-mass IPN parameterization, [90] for a point-mass 2PN parameterization, [91] for a quasi-Keplerian 3PN point-mass 
parameterization. [92, 93] for point-mass parameterizations under leading order spin-orbit coupling, [94] for using a 3PN point-mass 
parameterization including radiative dynamics for the phasing of gravitational waves, [95] for a post-equal-mass parameterization of 
a binary at 3PN point-mass level under leading order spin-orbit couphng, and finally [96] for a pai'ameterization up to 2.5PN with 
orbital angular momentum aligned spins. [97] will incorporate the linear-in-spin Hamiltonians given in the present article into the 
orbital elements for orbital momentum aligned spins. 
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series. Such a resummation was successfully implemented into the effective-one-body approach, which an- 
alytically provides complete gravitational waveforms for binary inspiral that are in good agreement with 
numerical relativity. This succeeded so far for point-masses [127, 128] and for non-precessing spins [129] 
by calibrations to full numerical simulations, but more work is needed for precessing spins [130, 131]. Here 
the Hamiltonians derived in the present paper should be very useful, and the spin-orbit one was indeed al- 
ready incorporated into the effective-one-body approach [132, 133]. See also [134] for a very complete 
overview of the hterature on the effective-one-body approach. Alternative ways of resumming the post- 
Newtonian series by Pade approximants are possible, which is most interesting for certain gauge invariant 
quantities [135, 136]. Within the overlap region of post-Newtonian approximation and numerical relativity 
in which the gravitational field is not too strong and the number of orbits can be handled by numerical sim- 
ulations the results of both approaches can be compared. The mentioned resummation methods can make 
these approximate results competitive to numerical relativity^ also in the late inspiral phase [127, 128]. 
Last but not least, a powerful interface between self-force calculations and the ADM Hamiltonians derives 
from the redshift observable and the first law of binary dynamics [150], which was extended to include 
spin recently [151]. 

For all computations we used xTensor [152], a free package for Mathematica [153], especially be- 
cause of its fast index canonicalizer based on the package xPerm [154]. We also used the package xPert 
[155], which is part of xTensor, for performing the perturbative part of our calculations. Furthermore 
we wrote several MATHEMATICA packages ourselves for the various steps of the computation including 
evaluating integrals. 

The article is organized as follows. Sect. 2 shows how to split the spacetime in a 1) manner, derives 
the Hamilton constraint and the momentum constraint from the Einstein-Hilbert action, and shows the def- 
inition of the ADM Hamiltonian. In Sect. 3 the constraint equations are expanded using a post-Newtonian 
power counting scheme and several integrations by parts simplify subsequent integrations. Afterwards a 
transition from the ADM Hamiltonian to a Routhian via a Legendre transform is performed in Sect. 4. The 
integrands are composed of various fields and their sources. These are explained in detail in Sect. 5. From 
the sources one can obtain the solution of the constraint equations and wave equations in Sect. 6 and is 
ready to integrate. Due to its importance a detailed explanation of the ultraviolet analysis is also provided 
there. In Sect. 7 the resulting next-to-next-to-leading order Hamiltonians are given. The kinematic con- 
sistency check of the Hamiltonians, namely the global post-Newtonian approximate Poincare algebra, is 
discussed in Sect. 8. Also the center-of-mass vectors are uniquely determined from an ansatz in the same 
section. Another check is provided in Sect. 9. There the Hamiltonians from Sect. 7 are compared with the 
Hamiltonian of a test-spin in a stationary exterior gravitational field. After that in Sect. 10 the conclusions 
are presented and in the Appendicesmore details on some of the calculation procedures used in the former 
sections are provided. 

The spacetime has d spatial dimensions, 1 time dimension, and metric signature d— 1. In this article a 
restriction to d = 3 is explicitly written, if this is not the case d is always generic. All calculations at the 
level of the field equations are performed in arbitrary dimensions, because of the necessary ultraviolet(UV) 
analysis concerning the short-range decay of fields around the sources. Vectors are written in boldface 
and their components are denoted by Latin indices from the middle of the alphabet. The scalar product 
between two vectors a and b is denoted by (ab) = (a • b). Our units are such that c = 1. There is 
no special convention for Newton's gravitational constant G^'^\ (Notice that G*^"^^ is the d-dimensional 
coupling strength of the gravitational field. It has the same numerical value as G in d = 3 dimensions, 
but different units.) In the results and the expressions for the sources, Pa denotes the canonical linear 
momentum of the ath object, the canonical conjugate position of the object, nia the mass of the object, 
Sa and Sa(i)(j) the spin vector and the spin tensor of the object, rat = |zq — Zfcl the relative distance 
between two objects, and Uab = (za — ^b)/rab the direction vector pointing from object b to object a. Also 
important are the distance between source a and field point, = |x — Za|, the unit vector pointing from 
source a to the field point Ua = {x — Za)/ra, and the circumference of the triangle of source a, source b 
and the field point x, given by Sab — Ta + rb + Tab- In the binary case the object labels a, b take only the 
values 1 and 2. The round brackets around the indices of the canonical spin tensor Sa indicate that its 

^See e.g. [137-139] for reviews, [140] for the first simulation of binaiy neutron stars, and e.g. [ 141-143] for very recent studies 
about them. The first simulations of coalescences of binary black holes were performed in [144, 145]. Furthermore see e.g. [146-148] 
for numerical simulations of a system of more than two black holes and the recent publication [149] about its chaotic behavior 
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components are given in a local Lorentz basis, which is essential for the canonical formalism, see [75, 76]. 



2 ADM Canonical Formalism 

In the following we introduce the ADM canonical formalism [22, 156, 157]. We pay special attention to 
the dependence on the spatial dimensions d, as dimensional regularization is important for the consistency 
of the post-Newtonian approximation when point-like sources are utilized [31]. The ADM formalism is 
extended from non-spinning to spinning point-like sources valid to linear order in spin here, see [75, 76] 
for the case d = 3. 



2.1 Setting the Canonical Formalism 

Canonical methods usually require one to single out a time coordinate, thus splitting spacetime into space 
and time. The geometrically favored way for such a splitting gives rise to the line element 

ds^ = -N^df + 7y {N^dt + dx')iN^dt + dx^ ) , (2) 

which corresponds to a (d + 1) -dimensional metric tensor field given by 

where N is the lapse function, A^' the shift vector, Ni — jijN^ and 7^ the metric of the spatial slices, 
7*'' being its inverse [158-160]. Notice that rewriting the metric tensor field in terms of lapse, shift, and 
the d-dimensional metric 7^ corresponds only to another representation of the metric, since together they 
have the same number of degrees of freedom. On the one hand a symmetric rank two tensor field like 5^1, 
in d + 1 dimensions has {d + l){d + 2)/2 independent entries and on the other hand, a symmetric rank 
two tensor in d dimensions like 7^ has d{d + l)/2 degrees of freedom, the vector N"^ has d independent 
entries, and the scalar N represents one degree of freedom. Obviously the degrees of freedom match. 
The action of the gravitational field is given by the usual Einstein-Hilbert action, namely 

Wfield - J d'^+lx Afield , /:field = ^^^^^'^'^^ ' (4) 

where g = det(5py) and (''+i)r is the {d + 1) -dimensional Ricci scalar, which can be split in a (d + 1) 
manner resulting in 



Add = Y^^N^ [(^)R + K,,K^^ - h^,K^'? 



(td). (5) 



Here '■''■'R is the d-dimensional spatial Ricci scalar, Kij is the extrinsic curvature, 7 = det{-/ij), and (td) 
denotes a total divergence which we ignore for now (it will be discussed in Sect. 2.4). We define 

fir 

n'^ = 167rG('') - V^i^'l"' - f''l'')Kki , (6) 

07^,0 

where 2NKij = — 7ij,o + "^^{i-j) was used (which deviates from the convention for Kij used in [159]), a 
comma denotes a partial derivative, and a semicolon denotes the d-dimensional covariant derivative. The 
inversion reads 

and the dimension finally enters the calculation explicitly. In order to put the field equations in the form of 
Hamilton's canonical equations, we perform a Legendre transform, 

CmA - T7r^^''l^lfi - + N^nf'^ + (td) , (8) 
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where 



n 



field 



field 



ik jl 



7 h^j^' 



(9a) 
(9b) 



Variation of the action with respect to and tt*^ leads to d{d + 1) first order evolution equations for 
these variables. The d-metric 7^^ and the field momentum Ye^^^^jW"^ canonically conjugate variables 
in the vacuum case. However, the action does not contain time derivatives of N and iV% so a variation with 
respect to these variables leads to constraint-type equations (i.e. containing no time derivatives). Variation 
of N gives only one equation called the Hamilton constraint, while variation of gives d equations called 
momentum constraint. We have arrived at a canonical formalism without any restriction on the coordinates, 
or gauge, but in the presence of constraints. In the following sections we will eliminate these constraints 
by simultaneously fixing the gauge. 

In order to couple the ADM formalism to point-like objects possessing masses nia and (d + 1)- 
dimensional spin tensors Sa = —Savfi, it is best to start with a matter action of the form 



E 



(10) 



where Tq is a worldline parameter, the [d + l)-momentum, — dz^/dTa, zj^ the position, 51^"^ = 
— rij^'' the angular velocity tensor, and Aq a Lagrange multiplier belonging to the mass-shell constraint 
9^'^PaiiPau + = (all for the a-th object), see [75, 76] for details. Notice that the angular velocity 
tensor is built from a Lorentz matrix encoding the orientation of a body-fixed frame and its co variant 
derivative. This means one has to handle a derivative coupling of the metric because of Christoffel symbols 
or Ricci rotation coefficients appearing in the co variant derivative of the Lorentz matrix. Further constraints 
must be fulfilled for the mentioned Lorentz matrix and the spin, the latter reads S^'^pa ^ = 0. It is important 
that this action is valid to linear order in spin and is now considered for a generic spatial dimension d. The 
equations of motion following from this matter action are the well-known Mathisson-Papapetrou equations 
[161, 162], see also [163, 164], 



dTa 



(11) 



where ^'^'^^^Rfj.ppa is the {d + 1) -dimensional Riemann tensor, and the source of the gravitational field 
equations is given by Tulczyjew's singular energy-momentum tensor density [163] 



(12) 



where (5(d+i)a = (5(a;^ — z^). For a review of spin in relativity see e.g. [165, 166]. Further details on a 
corresponding action principle can be found in [ 167-17 1 ] . 

Next the matter constraints are eliminated from the action with suitable gauge choices, e.g. Ta = t 
where t is the coordinate time and the matter variables are transformed to (reduced) canonical variables 
denoted by a hat, e.g. or pai- This is completely analogous to [75, 76] and will therefore not be repeated 
here. But the following facts are important for the present article. The spatial dimension d is not entering 
the just mentioned calculations explicitly, whereas it appears in the gravitational part of the action. Further, 
the matter action becomes linear in lapse and shift, so the gravitational constraints following from their 
variation now contain matter source terms 7^™^"'='' and 



n 



field 



0, n 



field 



n" 



0. 



(13) 



Finally, due to the coupling of the spin to derivatives of the metric, time derivatives of the gravitational 
field are present in the matter part. This necessitates a matter contribution Tr^atter to the canonical field 
momentum, which now reads 



(14) 
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Explicit expressions for T^"""":'^ 'H"^"'''', and Tr^atter provided in Sect. 5 and have the same form as for 
d = 3in [75,76]. 

The next goal is the elimination of the field constraints within the so called ADM transverse-traceless 
gauge (ADMTT). The corresponding gauge conditions read 

7ij J - ^7ij,i = , (15a) 
7r" = 0. (15b) 
We proceed to work out field decompositions suitable for an elimination of the field constraints. 



2.2 Metric Decomposition 

The differential gauge condition (15a) for the metric is solved by 

= ^S,, + hj] , (16) 

which can be seen from a transverse-traceless (TT) decomposition of 7^ . The first part is the conformally 
flat part of the metric and the last part can be interpreted in the far-zone as the radiation field, which is 
transverse-traceless, i.e. hH = = ^I/j - Due to the requirement of maximal simple curvature density 
expression [172], one can determine the conformal part of the metric in a form appropriate for a post- 
Newtonian expansion. Consider the static (i.e. momentum independent) part of the Hamilton constraint 
using (9a) and (13) setting tt*^ = 0, 

V7 ^''^R = IGTrG^'^^H™"" . (17) 

If one inserts the truncation 7^ = ^'(5^ with the ansatz ^I^ = ?/;^ into the static Hamilton constraint, one 
obtains 

-\p{d - l)Vj-2+5/3(rf-2) (4^AV' + (-4 + I3{d - 2)){'ilj^if) = 167rG('^)-H™"'='' , (18) 

where A = didi and di denotes a partial coordinate derivative. Demanding that the nonlinear gradient term 
(■^^i)^ disappears (yielding a Poisson-type equation) or the ip in front has a vanishing exponent, both leads 
to 

and thus gives the most simple expression for the curvature density '^''■'R. Using this solution, the static 
Hamilton constraint reduces to 

-lilzil^Ai/- = 167rG('')-H™"" . (20) 

Now one can further set ijj = \ + acj) and demand that the static part of the Hamilton constraint linear in (f> 
reduces to a Poisson-type equation without any d-dependence. This leads to 



d~2 

and 



(21) 



1 + f', ^, (/> ) A4> = WttG^'^^H""''''' . (22) 



4(d-l) 

Finally, the optimized metric decomposition reads 

d-2 



7«.-(l + i(rf3I)^) ^v+hll, (23) 
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see, e.g., [31]. One can also borrow the conformal factor from the d-dimensional isotropic Schwarzschild 
metric given in, e.g., [173, Eq. (22)]. For convenience we introduce 



in order to absorb certain dimension-depending factors in subsequent calculations. 



(24) 



2.3 IVIomeiitum Decomposition 

For convenience we introduce some differential-integral operators called 2?-operators, given by 

T^'^ = + (« - 1)5.9, A-1 ^ , (25) 
where A^^ is the inverse Laplacian with usual boundary conditions. These have the following nice prop- 



erties 



n 



and further 



Tf^ - ()■ ■ 
%2?j/' = d + a - 1 , 
d.V'l = ad, . 

Note that for a = one gets the transverse projector 

\ij — JJ'^J 



(26a) 

(26b) 
(26c) 



(27a) 
(27b) 
(27c) 



(28) 



which has vanishing divergence and is not invertible. From the multiplication property of the P-operators 
the projection property — follows immediately. Also a longitudinal projector can be con- 

structed, namely 



(29) 



Obviously the longitudinal projector also fulfills the projector condition = According to [76, 77] 

one can decompose the field momentum as 



The different parts are given by 



Tjj — 



J_k(ij_j)e 



1 



ij I ki 



d- 1 



^ij ^ I ij r ki eTrfe£ kl 



„kl . x^klM 



1 



d-1 



^ke gLJ k£^ke 



(30) 

(31a) 
(31b) 
(31c) 



Obviously tTjj has no divergence and is trace free (it is transverse-traceless), tt*^ is trace free and its 
divergence contains the divergence of the whole field momentum, and tt is divergence free and its trace 



8 



is the trace of the whole field momentum. This decomposition is essentially the usual transverse-traceless- 
decomposition of a symmetric rank two tensor field, but rearranged in a way more suitable for the present 
computations. For example, tt*^ is fixed by the gauge condition (15b) as 

1 I kk (^2) 

"matter ^ w^-l 



d-1 



which follows from its definition and the trace of Eq. (14). Furthermore one can decompose tt*-' into two 
different vector potentials, tt* and V^. The decompositions read 

^ 1^ti^%, (33) 



= V^+V^- ^6,,Vt . (34) 

Of course tt* and are interrelated, 

r = V%_, , = , F , (35) 

and it holds 

= A-^n]j . (36) 

These two vector potentials have different advantages. From it is possible to compute tt*-' without 
any integration. On the other hand tt^ has a simpler structure and can more easily be obtained from the 
momentum constraint using (36) [cf. (113), (1 14), and (115)]. The transverse-traceless parts of metric 7^ , 
hJJ and of the canonical field momentum tt*^ , tt^ are denoted as transverse-traceless degrees of freedom 
or propagating field degrees of freedom in the following. The latter denotation will become clear in the 
next section. 



2.4 ADM Hamiltonian 

We are now principally able to solve the d + 1 constraint equations (13) for the d + 1 variables and 
7r\ Though this involves solving a system of nonlinear partial differential equations, it can be tackled 
analytically within the post-Newtonian approximation up to a certain order. Notice that tt*-' is fixed by (32) 
and that tt^ can be replaced by 

"TT ~ "TT "ij "matter- \-" > 

The independent variables are thus the reduced canonical field variables hlj and irl^ with Poisson brackets 

{h^J(^), TT^^ (x')} = 167rG('^)<5TTfc.^(^ _ (38) 

as well as reduced canonical matter variables which enter via the matter parts of the constraint equations 
(13) and are discussed in more detail in Sect. 5. 

It was shown in [22, 156, 157] using different methods that the fully reduced Hamiltonian after gauge 
fixing is given by the ADM energy 

-Badm = ^Q^Q^d) / ^'^"^^Alr],] - In A , (39) 

where f d'^^^Si denotes an integral over the asymptotic boundary of a spatial hypersurfaces at fixed time. 
(The identical expression for the energy follows from the Landau-Lifshitz superpotential which is related 
with the well-known Landau-Lifshitz stress-energy-momentum pseudotensor of the gravitational field, see 
[174].) More precisely, the ADM energy Eadm turns into the ADM Hamiltonian Hadm if it is expressed in 
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terms of the mentioned reduced canonical variables. Inserting the metric decomposition (23), the surface 
integral can be transformed into a volume integral, 

where the asymptotic behavior of certain field components was used, see e.g. [ 157]. It was argued by Regge 
and Teitelboim [157] that one must modify the total divergence in (5) such that it leads to the surface integral 
expression for the ADM energy Eadm, otherwise the variational principle is not well posed. Indeed, for the 
same reason one should add the (regularized) York-Gibbons-Hawking, or "trace K," surface term [175- 
178] already to the Einstein-Hilbert action. The ADM energy then directly arises from a surface term 
contained in the complete action, see [179, 180] and also [159]. Notice that the Einstein equations can be 
followed from a variation of the action by disregarding all surface terms, but this is in general not allowed 
if the variation of the action has prescribed nontrivial behavior at the boundary. 

The ADM Hamiltonian still depends on the reduced canonical field variables hJJ and tt^j. In Sect. 4 
also these remaining field variables will be eliminated through a Routhian approach to arrive at a conser- 
vative matter-only Hamiltonian. 



3 Expansion of the Constraints and Integrations by Parts 

As already stated in Sect. 1 we will use the post-Newtonian approximation throughout this article. Further- 
more as mentioned above we utilize the spin-extended ADM formalism. So first of all we have to solve 
the constraint equations order by order This requires to expand them in powers of the post-Newtonian 
approximation parameter (1) with the decompositions (23) and (30) (which are adapted to the ADMTT 
gauge) inserted. This is the task performed in the following subsection. 

3.1 Order Counting 

The field and source expansions starting at their leading order are given by 

(41a) 
(41b) 
(41c) 
(41d) 

where the subscript in round brackets denotes the order. A similar order counting is also valid for all 
derived fields, like vector potentials. At a later stage of the calculation we also need to expand hJJ and ir^, 

hJJ - hf^y^, + hjly^j + . . . , (41e) 

%T = '^(5)TT + • • • • (41f) 






= 0(2) + 


0(4) + 0(6) 


+ 0(8) ^- 


- 0(10) + • 






- TT*-'' + 

- ^(Z) + 


^(5) +^(7) 


+ ... , 






^matter 


-27 matter 
- rt(2) 


1 nj matter i 
+ "-(4) + 


nj matter 
"(6) - 


1 'T/ matter 
^ ^(8) " 


1 -T/ matter 
^ ^(10) 


^matter 


•T/ matter 
- «i(3) 


1 nj matter i 

+ rLi(^b) + 


1/ matter 
«*(7) - 


h... , 





At the mentioned stage we also need an order counting for the deviations between field momentum and 
canonical field momentum due to Tr^^tter ™d the tracelessness violation of the field momentum tt . As one 
can see in Sect. 5.3 TTmaner starts at ©(c^^) and thus cannot contribute to the expansion of the Hamilton and 
momentum constraint later, while Tr^atter starts at 0{c^^). 

In anticipation of the source expressions in Sect. 5 we will introduce the matter variables used there 
to discuss the order counting of the field variables, because all fields depend on the source expression and 
thus the order counting of the matter variables. The mass ttIq, canonical matter momentum Pa, and spin 
variables Sa are formally counted as nia ~ 0{c~^), Pa ^ 0{c^^), and So ^ 0{c^^) for dimensional 
reasons only (remember that for maximal spins one would have ^ 0{c^^) instead, see e.g. [77, 
Appendix A]). This counting comes from the fact that after setting c = G'^'^^ = 1 we require all quantities 
to be in units of length. Let us introduce symbols with a bar below them being the quantities in SI units 
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and the other symbols the quantities in units of length, then it holds nia = ^^^ELa for 'he mass, t = ct 
for the time, Pa — ^=73- for the linear momentum, and similarly for the spin variables. Although we 
mentioned that G^''' has different units in d 7^ 3 than in d = 3, we treat their order as in d = 3 for 
simplicity. So the order counting comes from the c powers inserted to reconstruct the SI units. It should 
be noted that these counting rules wiU in general not give correct absolute orders in c if the SI units of 
the final expression are not taken into account. However, relative orders are always meaningful, which is 
all that is relevant for perturbative expansions. Further notice that different counting rules are obtained if 
one assumes that all quantities are expressed in terms of mass units instead of length units when setting 
c = G^*^^ = 1, which is also often used in the literature. 



3.2 Hamilton Constraint 



The Hamilton constraint 
1 



(42) 



167r^G('*) 

has to be expanded in a post-Newtonian manner After inserting the metric and momentum decomposition, 
(23) and (30), and all field expansions (41) it decomposes into several Poisson equations for the post- 
Newtonian potential, namely 

1 



167rG('*) 
1 

167rG('i) 
1 



1 

167rG('^) 



A(/'(2) 


nj matter 
- «(2) > 








A(/)(4) 


nj matter 

- n.(4) - ( 


1 o; matter 
P{2)Tt{2) ' 






A(/'(6) 


nj matter 
- «(6) - ' 


/'(2)^r4r+ 


(-0(4) + 


0(2) )rt(2) 



1 



167rG('^) 



+ 4(0(2) ) , 
matter 



(43a) 
(43b) 

(43c) 



(8) 
3n 



2)^(6) 



0(2) )H 

4 



matter 
(2) 



0(4) + 0(2) + (-0(6) + 20(2)0(4) 

3d - 10 



16 - - XT - XT 

^(0(2)0(2)/^^),. + 4(0(4)/i]7),., 

(43d) 
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-^A0(io) = H^-f - 4>i2)n^sT + (--^(4) + 0(2)')^r6r + (--^(6) + 20(2)0(4) 



0(2) jrt( 



■•(4) 



+ 20(2)0(6) + 0(4) - 30(4)0(2) + 0(2) 



(hi 



TT \2 \ nj matter 



1 



(2) 



3d - 10 



167rG('^) V rf - 2 
2- 



(^5)) - 2V3)Y7) - (%t) - 2 ur~W2 (Vs)) '^(2) 



TTr-ife ~jk 
y '^(3)'^(3) 



+ 



id-2r 

^ _ 0(2)^(4)^^- + 2 ^(3)^(5)^(2) + 2 ^3)%T0(2) 

,^^L hTJ$(2^ :4>C2\ _.0«^ - \hJJ j^h^l - \'K-^{h]lkf'^{2) 



Xd -"2)2 "^J" '^(2), ,9(2) ^^-9(2) 

+ (td). 



(43e) 



By virtue of (40) the post-Newtonian Hamiltonians follow from an integration of the right hand sides of 
these equations. The last equation leads to the formal 3PN ADM Hamiltonian. 
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3.3 Momentum Constraint 



Also the momentum constraint has to be expanded in a post-Newtonian manner. First of all one has to 
write the co variant divergence in a more explicit form. For convenience it is also useful to write as much 
terms as possible in terms of divergences of a traceless symmetric tensor, see [77]. Notice that we also did 
not remove the trace parts of the field momentum in 



1 



1 



j4/(d-2) 
1 



d-2 



1 



(44) 



as the tracelessness condition may be violated for the non-canonical field momentum due to the spin. The 
post-Newtonian expansion is necessary for the later integrations by parts and to obtain the necessary field 
solutions for the momentum type fields. The expansion reads 



^(3) J 



^(5) J 



+ vL h 



TT 

(3) I '^jk^t 



2(rf-6)., 



, TT _ , TT 

"■ifcj' ij,k 



(45a) 
(45b) 

(45c) 



Notice that we did not insert the expansion for hfj and vr.^^ since this is only possible after their evolution 
were obtained and solved order by order later on. 



3.4 Integration by Parts 

Due to the complicated structure of — A(/)(io)/(167rG'^''') in particular the appearance of 0(8), 0(g) and 



'(5)' "(7) 



it is necessary to simplify the integral over the right hand side of (43e). Some of the mentioned 
fields are not even known in d = 3 dimensions. The best way to remove them is to integrate by parts certain 
terms and afterwards use lower order Hamilton and momentum constraints. 

In the used dimensional regularization one may always neglect boundary terms if the integrands are not 
UV- and IR-divergent simultaneously. These more subtle terms occur at 4PN point-mass calculations for 
the first time. In the present calculation we always neglected boundary terms in the integrations by parts. 

The parts of the Hamilton constraint coming from the expanded Ricci scalar in the conformal approxi- 
mation have always a structure where a power of the post-Newtonian potential is coupled to a matter source 
of the Hamilton constraint. These terms can be simplified by inserting the lower order Hamilton constraint 
for the source and shift the emerging Laplacian to the coupled post-Newtonian potential via integrating by 
parts twice. This procedure can be used to eliminate 0(8) and 0(g) and the appropriate calculations are 
given by 



o_/ matter 

«(2) 



2 7 . 7 4^^, matter ^ J '^'^ ~ I ^ f ~ ij \2 

^ -^'2) (7^(3)) 



d-2 



- 20(2) 0(4) + 0(2) )n1^) " - ^Q^Qid) I 



16 



+ 7^(0(2)'^(2),»/l|/')j'?i'(2) - 4(0(4) J 0(2) - ^0(2)A(/lJ/)2 
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I ^ / matter 
-'P(6)rt(4) 



+ 2</'(2)(/iifeMfc),^j| + (td), 

7 ^ /matter > J J n / matter \ ( ~L 7 7 ^7 \n /matti 

-0(42) rt(6) + 0(2) 0(4 2) rt(4) + (0(4) 0(4 2) " 0(2) 0(4 2) ) rt(2) 

Y^^{0(4 2)(^;'3))' - 40(4 2) (0(2) /i^^),., I + (td) , 



30(6)0(2)H^2f 



-3H^;;f ■ + 30(2)H^^f ' - (-30(4) + 30(2)OHr2T 



1 



3(^;^3J2- 120(2), ../jTT 



0(4 0) + (td) . 



(46) 



(47) 



(48) 



The fields 0(4 2) and 0(4 q) are the momentum dependent and the momentum independent part of 0(4) given 

= -16^G('^)4^A-iH™"" and 0(4 2) = -\^t^G^^^ It is possible 



by 0(42)- .^-.^ 4^31)' 

to simplify (46) once more using 



1 



167rGW 



(2)A(/.]7)' 



d-2 



8{d- 



(49) 



In the Hamilton constraint there are also several terms coupling different orders of the longitudinal 
field momentum. These terms can also be simplified by inserting the decomposition (34), removing the 
derivatives from the vector potential via an integration by parts and inserting the lower order momentum 
constraints (45) into the divergencies of the longitudinal field momentum. This procedure is used in order 
to ehminate tt^^^ and tt^^-!^^ . The respective integrations by parts are given by 



^ ( ')~i3 ~ij \ — oxri o/matte 



1 



167rG('') 



-2h 



+ rf7-^^3)('^TT + *(5))0(2) - 4 

+ ^(^^3))%)K(^d). 

Last but not least the (7r(5))^ integration by parts is given by 



"[-2^(3)/(3) + (^(3)-r3)),^ 

d — 6 ,9 T 2 



(d-2)2^ (3) 



(*f3l)''^(2) 



167rG(^) V ^ (5)^ ) ~ ^(5) ^'(5) 



(td). 



(50) 



(51) 



16^GW\rf7:2"(3)"(5)<^(2) 

Although it is good to express field integrals in terms of sources and fields, the V^^^ potential is still very 

complicated. Thus we try to express (7r(5j)^ in yet another way. From the transverse-traceless projection 
of an arbitrary second rank tensor field , namely 



A,... 



Ay - A 



(52) 



(53) 



one can get the transverse-traceless projection (0(2)7i^?3\ is traceless) 



^IJ'' (^(.2)^(3)) = 0(2)^(3) + ^-J^ {^(5) + ^5)l} 

via the momentum constraint (45). Here 

1 



'-IJ -t 

'^{5)1 - ^(5)l,i 



+ ^(5)1, 



Ti^j ~k 



(54) 



(55) 
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where 



— 1 17 matter 
»(5) ■ 



(56) 



Notice that tt^^^j^ is defined with a different sign than the usual vector potentials used throughout this 
article. Now we can express yr^g^ in terms of tt^s)]^, 0(2) '''(3) and the transverse-traceless projection (54). 
From these considerations it follows that (7r|g^)^ is given by 



1 



167rGW 



16 



167rG('^) \{d- 2)2 



d-2 



(57) 



Some of the results above were partially checked by comparison with [77]. 

3.5 The formal 3PN ADM Hamiltonian 

Performing the mentioned integration by parts leads to the following formal 3PN ADM Hamiltonian which 
can also be compared to [77], 



H 



3PN 



'^O-o'r ~ 2(?!)(2)'H™"''' + (-5(4) - 40(4) + 20(2) ^)'H™)"" + (60(2)0(4) + 5(4)0(2) 



2 



20(2)')^;^)'" + ( 40(4)^ + 5(4)0(4) - 80(4)0(2)^ - 5(4)0(2)^ + 20(2)^ 



4(d-l) 



/ L 1 1 \ ^ \ nj matter i o T/* -IV matter 

y(/''y-J j«(2) +^^(3)«»(7) - i67rG('i) 



-(-r5)i)' + 2^^(4)(-r3))' 



+ %)(-r3))' - '^'\d% '^ (^3))^^(2)' - (4^)^ + 80(2)V^;^3)4^ 



2W 



2^(^)/g) + (^(;^)-r3))-^ + '^(3)-(3) - 



k ~ij 



^ik ~jk 



.2d -3 



■0(2) ,,0(4), 



3d -4 - 

2d -3 
{d^ 



(2)'P(2).,'P(2), 



25, 



16 



(4),jn2)j 
2 2 



rf- 1 



8^^0(2)^(5)1^3) 



^fcr^(0(2)-f^))j +^(^;^)'0(2) 



(58) 



We changed all occurrences of 0(4 2) to 5(4) = —20(42) to gain a result which is comparable to [31]. Note 
that, since we did not expand hJJ and tt^, we also need some formal 2PN terms which may contribute to 
the 3PN kinetic energy or the 3PN interaction Hamiltonian after expanding hJJ, 



nj matter 
''•(8)TT " 



1 



167rGW 



4(d-l) 



+ ^i:^0(2),^0(2)>"-^(/^"fc)' 



(59) 



Now we can split up the Hamiltonian into a kinetic part for the reduced canonical field variables [hJJ and 
TTjj, after inserting (37)], an interaction part between these canonical fields and constraint fields, and a part 
independent of the canonical fields, i.e.. 



ADM — -^ADM -"ADM 



int 



H 



non-TT 



1 



167rG('^) 



(60) 



with 



1 



^™ 167rG('^) 



d-x\ (i?(4)„- + B(,),,)hjj - ^^(/,TT)2^n.atter ^ ^^^.^ 
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2 



0(2) (/l. 



rrnon-TT rr 



TTint 

-"ADM 



1 



and 



jd^ -T/ matter 
^^TT / " "-(8) 



4(rf-l) ;, 7 

'P(2),»'?(2),J , 



(61) 
(62) 

(63) 



B 



(6)»J 



167rG('')- 



u 



(10) 



O ^ -IV matter , o T / -IV matter 



^2d- 3^ 
d- 2 ' 



3d -4 



(2), . -^(4) + 4-^^^3^0(2)0(2) ,,0(2), 



- 2%) ^.0(2)^ 



^"matter 



80(2)^(3) 



(64) 
(65) 



^ADM consists of all Hamiltonians starting from the rest mass contribution up to Hspn- The B(4)ij part in 
-^ADM comes from the formal 2PN Hamiltonian shown in (59). U^^^^ can be obtained by removing all 
and parts in (59) and (58) respectively. Notice that the relevant source terms in the expressions for 
B(4)ij and B{%)ij are at most linear in li^ . This allowed us to single out these contributions by a functional 
derivative. 

Now we arrived at a point mentioned in Sect. 2.4, namely where we are able to eliminate the constraint 
fields using lower order Hamilton and momentum constraints, but where the dynamical field degrees of 
freedom are still in the Hamiltonian. The elimination of these and further simplifications of the calculation 
process are subject of the following section. 



4 Routhian and Application of Wave Equation 

To obtain a fully reduced matter only Hamiltonian, we have to remove the dynamical degrees of free- 
dom /iJJ and TTj^ by solving the appropriate equations of motion and inserting their solutions into Hadm- 
However there are some subtleties in this procedure which will be discussed in Sect. 4.2. 
From Hamilton's equations 



1 



167rG('') 
1 



'^^IJ _ cTTij ADM 



at 

^ ^ " A/lTT 



167rG('') dt 

and the split of the ADM Hamiltonian (60) one gets the appropriate wave equations 



(66) 
(67) 



ahJJ = i6^G('^)j,, 



TTij 



s: zijint 
o '^-"ADM 



9 (^-f^ADM 



dt Sttjj 



'''tt ~ 



r ijint 

167rG('^)<5^^'^'^:^ 



(68) 
(69) 



for the dynamical degrees of freedom of the gravitational field (here □ = A — c^^d^; remember that we 
use a {d + 1) -metric with signature d — 1). In order to receive the full wave equations one has to perform 
the variational derivative of the interaction Hamiltonian, 



ah 



2(i?(, 



B 



{6)ke) 



levrG^-^) 
d- 1 



nj matter j TT 
''•(2) "■kt 



d-2 



(0(; 



?,TT N 

)"'kl,m),rn ' 



d_ 

dt 



Cki 



(70) 
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hjl - - 

2 2 



hJJ--jZ''C,,. (71) 



4.1 Near-Zone Expansion 

The wave equation 

nh = f, (72) 

[which represents the components of (70)] for a source / has several solutions depending on the boundary 
conditions. In field theory mostly the retarded one is used. At the order considered in the present article 
the time symmetric (i.e. conservative) solution is sufficient.-' 

The appropriate solution derives from a near-zone expansion, which is a formal expansion in c^^. 
(see e.g. [181-184].) More precisely, the near-zone expansion is a series expansion in the small quantity 
r/{ct) <C 1, which means that the distance from the field point to the appropriate source point is small 
compared to the gravitational wavelength. So the near-zone expansion may be used if retardation effects 
are negligible. Consider the Feynman propagator for a massless particle (which corresponds to the Green's 
function of the wave equation) 

G^(x,t) = --lim /dfco-— . d'^k——^—. (73) 



2tt £->o J " {271^ J k^-k^-ie 

For the reader's convenience we reintroduce the powers of c in the following expressions. This means 
fco = w/c and a;° — ct. As it was argued above and in e.g. [185, 186] the condition r/ (ct) <C 1 corresponds 
to fco/fc ^ 1 which gives rise to the so-called potential gravitons in contrast to the radiation gravitons for 
which fcp « k^. For the radiation gravitons the ie term becomes important, but for the potential ones it 
could be neglected and the Fourier amplitude in the propagator is expandable in fco / k, namely 

^ ^ 2ttcJ ^ c2" {2nY J fc2(n+l) 

This means that a near-zone expansion of the (time symmetric) Feynman propagator cannot contain any 
contributions leading to radiative losses. The d-dimensional fc-space integral in (74) has to be performed 
by using 



(47r)d/2 r(a) V 4 

which can be obtained by dimensional regularization (the source of the wave lies at the origin of the 
coordinate system). This result can be reformulated in inverse Laplacians on a delta-type source by iterating 
(167) and (168) in the Appendix, 

1 f.dul^- r(2 + n-f)r(f-i-n 



{2TT)d 



fc2(n+i) rf2-|)r(|-i) ^ ^ 



For d ^ 2Z (i.e. no odd dimensional spacetime) one can simplify the Gamma functions further by using 
the identity r(l — z)r{z) = tt/ sin(7rz), z ^ Z which leads to 



d'k^^^-{~inA-r^'s, ill) 



{2itY J fc2("+i 



^^At the considered order ^J^-^^j '^^n be neglected since the linear terms are not time symmetric and the quadratic ternis arise 
only due to (hj^.^^^. j.)^ which is zero in d = 3 because '^IJjjj is only a function of t (there are no explicit terms in the source at 
c~^-order); ^T^-^^j is of too high order to appear at formal 3PN order, see [32, 59, 77]. We therefore neglected ^^IJ^jj ™d ^'(7)ij 
(41e) and hence there are no contributions to (58) or (59). 
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and therefore after performing the lo Fourier transform to 



NZ 

n=0 



1 d 



2n 



GH^,t) E ( ) (A-^)"+M(,+i) , (78) 



for the Feynman propagator This result could immediately be used to write down the near-zone expanded 
solution of the wave equation 



n=0 



1 O 



Notice that a near- zone expanded field in general does not converge at spatial infinity. 

Now we are able to derive the solutions for the transverse-traceless part of the metric at a certain post- 
Newtonian order in the near-zone. 



4.2 Routhian 

Before we can insert the wave equation (70) and its solution [see (79), and for a more explicit form see 
Sect. 6.2], we have to transform the ADM Hamiltonian into a Routhian, i.e. a Lagrangian in hlj and tt!^ 
and a Hamiltonian in the particle degrees of freedom, 

i?[/iTT, ATT] ^ ^^^^ _ _}_ J ^^^j^TT _ (gQ) 

This is necessary because otherwise the equation of motion, e.g., for z^, following from the Hamiltonian 
(for simplicity we omit the spin variables here) 

HADM{2a,Pa,hJJ{Zb,Pf,),TT!^{Zb,Pb)) , (81) 

using the Poisson brackets would be 

|fc _ dHADM _|_ SHadm ^^JJ _|_ SHadm d-k^ ^^^^ 
dpak ShJJ dpak dpak ' 

This equation is obviously wrong, because hJJ and ttjj are dynamical degrees of freedom and may not 
lead to additional terms in the equations of motion, also when their solutions are inserted. On the other 
hand using R{za, Pa, h]J{zi,, Pb), h]J{zb, Zb, Pb, Pb)) one has for the same equation of motion 




^a^^ + Tfrr ^ , (83) 



which has no additional terms coming from the chain rule, because they vanish due to hJJ fulfills the 
equations of motion in the appropriate approximation. Hence one can obtain the equation of motion in the 
usual way and one does not have to keep track of the field insertions. This is analogous to the construction 
of a Fokker action, see, e.g., [37] and references therein. 

A Fokker-like construction of a matter-only Lagrangian (or Routhian) can not account for dissipative 
effects, see [187] for a discussion. However, dissipative Hamiltonians in the ADM formalism can be 
constructed in the following way [32, 59]: The matter variables entering the solution of hJJ and tt^ 
are substituted by new ("primed") variables and thus the binary Hamiltonian now contains four types of 
canonical matter variables. This procedure prevents occurrence of wrong contributions in the equations of 
motion, too. The primed variables will be treated as explicitly time dependent and lead to an explicitly 
time dependent Hamiltonian. After calculating the equations of motion for canonical positions Zq and 
momenta p^, the primed variables are again identified with the old ones which makes another regularization 
procedure necessary [32]. Another approach to construct an action principle for dissipative systems was 
suggested in [187] recently. 
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4.3 Reduction of the Routhian using the Wave Equation 

The part of the Routhian labeled as TT contains iJ^DM' '^^e field kinetic part of the Hamiltonian and the part 
coming from the Legendre transform (and thus only terms coming from the transverse-traceless degrees of 
freedom), 



1 



167rG('i) 
2 - 



, TT x2 
ij-k) 



(4)iJ 
1 

4' 



i TT 
''ij.k 



Inserting OhJJ and tt^ from (70) and (71), or 







d-1 


TT,* 


"'I J "TT 



TT\2 nj matter 
«(2) - 



(84) 



hJJ{B(^i)i.j + B, 



(6)»J, 



1,TT_^^. . 



d-l 
-(td), 



/ 7 TT \ 2 -2 / matter 
J ^(2) 



d-2 '1 



(85) 



leads to 



1 



167rG('^) 
ld_ 
2dt 



TT 



tTT 



Ci. 



1 



— A 



47rG('') 

^ _ 1 "-"-^i > ^(2) 



LTT\2-2/matter 



Gfef 



d- 2 



(86) 



where the last part will appear in the matter part of the final Routhian, because there is no hJJ and no tt^ 
or h]J in Cij . Notice that we kept a total time derivative here. If we would drop it the Cij terms would 
not cancel in the next step. These terms are not impossible to handle, but it is advised to remove them to 
simplify the calculation. 



4.4 Insertion of the Near-Zone Wave Equation for Further Simplification 

Now we need to split up the first expression in the TT part of the Routhian (86). The near-zone expansion 
of the transverse-traceless part of the metric hJJ — hj^-^^ - + h^^^^ + . . . , which is explained in detail in 
Sect. 4.1, contributes only via /i(J)y and hj^-^^- at 3PN level. This expansion leads to 



where (ttd) denotes a total time derivative. Notice that there is a difference between □^(Jj and (□ft,JJ)(6) 
in the near-zone expansion as time derivatives raise the formal order of a field in contrast to spatial deriva- 
tives. This is a specific feature of the near-zone expansion, as in the far-zone time and space derivatives are 
equal in magnitude. These considerations lead to the difference in the following box operations. 



Ot "■(4) ij 



(88) 

(89) 



where dfhj^^ is of formal c * order and hence the total time derivative can be neglected. From this it 
follows that 

r TTniLTT 



"■(4) 



lTT a2LTT 
"■(4) Jj0'i"-(4) I] 



-(td) 



(90) 



such that 



y n^i7)(io) - 2/1(7)^- (□/i")(6) + hjjy^^dthjjy^j + (td) . 



(91) 
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where (□/iJJ)(e) is given by the c ^ part of Eq. (70). After another integration by parts this immediately 
leads to (we need only the leading order of Cij which is at c~^) 



jTT p ^fh^T: \2 ^^^^'^\i,TT ^2a7matter 



/1, 11 xznimdtt' 



1 



(92) 



Rearranging some of the terms into pure matter and pure trans verse-traceless parts, one obtains the full 
final 3PN Routhian, 



n 



(10)n 



TT - 24>i2)nf^^^,.r, + (-5(4) - 4^(4) + 2^(2) )H 



matter 
(6) 



+ (6<^(2)0(4) + 5(4)0(2) - 20(2) + I 40(4) + 5(4)0(4) " 80(4)0(2) " 5(4)0(2) 



+ ^0(2) I rt(2) + ^ ^(3) ^i(7)non-TT 



1 



167rG('*) 



-(-;^5)i)^+2^0(4)(*;^3))^ 



. s2 (3rf-4)(3d-2) , 27 
+ %)(Y3)) (^_2)2 (Va)) '^(2; 



2 , / 4(d- 1) ^TT»jVI 

^fc^ m2)^3). 



(3)^ 
2 



(93a) 



^3PN — 



1 



167rG('^) 



I'iTT \2 ^TtG^'') \2T/mattei- 



111 p ('f,ll ^-'^ " i^;.ll VZumitb 

"•(4)ij^(6)u - 4l"(4) jjJ ~ ^_ ;^ '>"'(4)ijJ «-(2) 



^^0(2)(/l")ij- fe)^ - ^('4^)^j'r(^5)matter " ^0(2) ^(3) ^('I)y 



(93b) 



.TT ^2 ;,TT _ A A. .^^i,%^h^^^ 

where Bi^Q^ij is given by (64). Note, that one does not need to calculate the hj^-^^ - field, which is not fully 

known in closed form. An explicit form of B{6)ij with all sources inserted is derived in the next section, 
cf. (107b). 



5 Sources 

The construction of the sources linear in spin follows along the lines of [75]. This requires the introduction 
of a local Lorentz frame"* as we want flat space Poisson brackets for the spin. In particular we also apply the 
Schwinger time gauge for the {d+ 1 ) -dimensional framefield which effectively reduces it to a d-dimensional 
spatial framefield e(i)j . 

During the following calculations we neglect the Tr^atter terms, which are far too high in their order, 
such that the modified source terms are given by [76, Eqs. (6.33) and (6.34)] 



t; matter _ „~ i; Pajl''^ ^kl „ (m) r 



A np^ ^ {npJHma-npJ P y + ^ "^y'"'' 



Pa.e ij kl s- \ 1 c 

——7-7 '7 e(r)je(s)k(>a\ \^a{r)(s) 
'"•a ii-Pa / J 

^matter ^ ^ _ A^^'e^^^J^Sa + \ f 7'"'e(,),e(,)fc<S, 



(94) 



"^Such frames were originally invented by Elie Cartan and named "repere mobile", which is French for "moving frame". In d = 3 
it is also called "triad" or "dreibein"; in d = 4 "tetrad" or "vierbein". For arbitrary integer d it is called "vielbein". 
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(95) 



where ttIq is the mass of the ath object, Pa its canonical momentum, Sa{r){s) its canonical spin, np^ = 
— y/ + 7*-'pa iPa j , and (5a = (5(x — Zq ) is the d-dimensional Dirac delta located atx — Za- Furthermore 

= f S^-^ (Is^,, + ) , (96) 

= -P!in!!l^ , (97) 

to linear order in spin. However the 4 -terms do not contribute to the expanded source expressions at the 
orders considered here. The matter position and matter momentum variables are canonical conjugate to 
each other, namely 

{zlpbj}^S,,Sab, (98) 

and the spin variables fulfill also canonical Poisson bracket relations, namely 

{Sa{i){]), Sa{k){i)} — 5ikSa{j)(l) — 5uSa(j){k) — SjkSa{i){e) + 5jtSa{t){k) , (99) 

where the canonical spin tensor is related to the canonical spin vector Sa via Sa = SijkSa (k) 

and Eijk is the Levi-Civita symbol. The appropriate Poisson brackets for the canonical spin vector are given 
by 

{Sa{i),Sa(j)} = £tjkSa(k) ■ (100) 

5.1 Framefield Expansion 

We choose to work within a symmetrical framefield gauge e(i)j = ei^j^i [188], so the dreibein can be 
written as a matrix square root of the metric, symbolically e(j)j — or more explicit 

e{i)ke{k)j = 7ii ■ (101) 

Notice that jij is positive definite and we require the same for e(j)j such that it is unique. The second 
relation, (101), can be inverted order by order, namely 

ADMIT e /im \ 

6(0) (i)fee(o) (fe)i = 7(0) i] 6(0) {i}] = Oij , (lOZa) 

ADMIT 1 2 7 , /imu\ 

6(2)(»)fce(o)(fe)j +e(o)(j)fe6(2)(fc)j =7(2) ij ^ 6(2) (,)j - -7(2) - ^0(2)Ojj , (102b) 



(103a) 
(103b) 



and at the end of the day 
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e(4)(.),- (^_2 
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2)2 




0(6) 


d- 
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e(6)(.),-(^_2 




2)2 



^)s^, + lhJJ, (103c) 



2(d-4)(d-3)- 3 



^(^)""3 (d-2)3 '^t^) /'^■-d^'^(^)'^'^ 

(103d) 

for the framefield perturbations. The antisymmetric part of the framefield (which is zero in this gauge) can 
be interpreted as rotational degrees of freedom in the choice of the local frame. Such a rotation does not 
change the length of the spins. Recall that an antisymmetric matrix is an infinitesimal generator of rotations 
and in d dimensions has ^ d{d — 1) independent entries. This is exactly the number of rotation planes in d 
dimensions. 



II 
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5.2 Constraint Sources 



After performing the expansion of the sources (94) and (95) in powers of c^^ and also expanding the 
metric-framefield relation (101), we are able to write down the appropriate post-Newtonian contributions 
to the source of the Hamilton constraint, namely 



'T/ matter 
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(104b) 
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128m7 " 4(d-2) ml 



<t){2)Sa 



2{d-2f ml 



4>{2) 5a 



2d{d + 2) pg 
3((i- 2)3 m„ 

2 pg 



0(2) (^a 



1 (pg)2 2(rf + 2) pg - 

'^(4)'^° + M 0^2 — '^(2)'^(4)Qa 



d — 2 ma 
3 (M)^ 



(d- 2)2 ma 



P^ 



0? — 2 to3 

i TT 



ff + 6 p2 - - - 

4(d- 2) „j5~P''*'^'^W0")'^(2),/« + (d-2)2;^^"'^°(*)(j')'^(2)0(2)j'5a 

2d(d + 2)23a^^ 2- 1 pg ^ c I A 

-'^a(»)(j)0(2) 0(2) - T^"!^ ,:ia (»)(j) 0(4) .^-^a 



(d-2)3 ma 

2(d + 2)pai 



d- 2m3' 

2 Pai 



(^i:^— ^a(»)(.)(0(2)0(4)),,'5a + ^ " (.)(,) 0(6) , /a 
PaiSa(j)(k)hJj^kSa- ■J—^^^Sa(j)(k)hJjj^(j)(2)5a 



Ami 



3 T^a z , TT 



d — 2 ma 



hikSa{j){k)4'i2),Sa 



d — 2 ma 

1 Pa, 



d — 2 ma 



'a{i){k)hJl(t)(2).jSa 



(td). 



In ^"o'" there is a 0(5) term left. But all occurrences of 0(6) can be cast into the form 
which we integrate by parts using (47). Then 0(g) disappears and gets substituted by 
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(td). 



(105) 



Furthermore we are also able to write down the sources for the momentum constraint in their full form, 
which are given by 
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(106a) 
(106b) 



A2 
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-^{hjkSa(i)(k) + hikSaU){k)) + -^;^Pak{PajSa{i)(k) + PaiSa {j)(k)) 



(106c) 



With the source expressions (104) including (105) and (106) from above, we may express B(4)y and B{e)ij 
in terms of the matter variables. They are given by 
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(107a) 



(107b) 



We expressed ^^(6)ij in ^ more convenient way now, since the tt vector potential has a much more simple 
structure than the V^g^ vector potential. 



5.3 Matter Correction to the Canonical Field Momentum 

Since we eliminated the transverse-traceless part of the canonical field momentum via using the relation 
between canonical field momentum and velocity of the hlj field (71), there are terms containing matter 
parts of the field momentum left in the Routhian (93b). These can be calculated from e.g. [77, Eq. (2.34)] 
where 7r^^j,tter is given by 



(108) 



From [77, Eqs. (3.33) and (3.34)] one gets the closed form expression 



^-3 ^ l..fc.,,-£ -n^aPa^kUSa,) 
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The part containing the antisymmetric Aa^^ was neglected, because it is of order c^^ (B^u starts at c"'') 
which is not necessary here. Here we only need tt^J to the order c~^. A power counting (see beginning of 
subsection 3.1) tells us that we only have to take the leading order approximation of the above expression, 
reading 

^5)a = ^ {pazSaik)U) + Pa 3 S a (k)(t)) ■ (HO) 

This expression has a vanishing trace, which makes obvious that we can neglect tt*^, Eq. (32), at the 
considered order. 



6 Field Solutions and Integration 

After obtaining expressions for the sources (104), ( 106) and the Routhian (93) which gives the Hamiltonian 
in the matter degrees of freedom after an integration, we need the fields to be inserted into the Routhian. 
These can be derived by solving the lower order constraint equations in case of the post-Newtonian potential 
and the non-propagating parts of the field momentum (see 6.1). For the propagating degrees of freedom 
the wave equation has to be solved (see 6.2). 



6.1 (i-dimensional Solutions of the Constraints 

rC^ — i) 

With K = d G'-'^\ the Hamilton constraint equations (43a), (43b), the momentum constraint equation 

(45a), and the various transformation formulas (33), (34), and (35) relating the longitudinal field momen- 
tum and its corresponding vector potentials, we find using the inverse Laplacians Usted in the AppendixA. 1 
that 
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(Ill) 

(112) 
(113) 
(114) 



(115) 



Remember that the momentum constraint can be solved for tt* with the help of (36). The more complicated 
fields like Tr^gj or (/)(g) were so far only found in c? = 3 dimensions [27]. Also the leading order of the 
transverse-traceless part of the metric is only partially known in d dimensions. We will discuss these issues 
in the following subsection. 



6.2 Solutions of the Wave Equation 

Consider now the wave equation (70) for hJJ. There hJJ is given in terms of a post-Newtonian approximate 
source Sij, namely 

^f^lj = S'ke''^i^(.i)ki + S(Q)kt) ■ (116) 
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By taking into account the near-zone expansion of hJJ, (79), one gets 



(4)y ~ "ki 

rTT 



(4)fe£ : 



^(6)y - ^ke'^^ ^ {S{6)ke 



(117) 
(118) 



for the leading order and next-to-leading order expressions of hJJ in the near-zone. Here the sources are 
given by 



S{e)tj ~ 2i?, 



167rG('^) 
{6)ij ; — —n 



(119) 
(120) 



where B(^4)ij and B^())ij are given by (107a) and (107b), and C(5)ij by (65) via (108) and (110). Notice 
that we removed the post-Newtonian order-counting parameter c in (117) and (118). Fortunately there 
is no need to evaluate (118) here. In fact the hJJ^^^ dependence of (120) renders the calculation almost 
impossible. 

Then the solution of the wave equation at leading order and linear in spin is given by 
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(121) 



where hj^^^y^^ is the momentum (and spin-) independent part of the transverse-traceless part of the metric 

which is generated by the TT-projection of A~^(0(2),i'^(2).j) ™d is only known in d = 3 see [27, Eq. 
(A20)], namely 

1 f ra+n 
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(122) 



Both solutions were also obtained by using the inverse Laplacians in the AppendixA.l. Obviously, most 
of the parts of hJ^-^^_^ are of the same type as hj^^-^^^ (see (118) and (120)). That is the reason why we 
eliminated it from the integrands. 



6.3 Distributional Contributions 

As long as the Riesz-kernel method is not used (where a Dirac delta is substituted by the so-called Riesz- 
kernel) one has to take care of delta parts when differentiating certain functions. Consider for example the 
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field 0(2) = 4:K J^a '^^'^ differentiate it two times. Using the ordinary derivative it would give 

dT'dYU{2)^o- (123) 

But as we already know from the constraint equations the second derivative of (j){2) should be 

^^^,(b(^2) = -167rG('^)i% ^ma(5, . (124) 

a 

Fortunately there is a result from the theory of distributions [27] which defines the so-called distributional 
derivative 

= Off + ^^^J^ £ d^^.-i n7.^^ ... .V (125) 

Here / is a positive homogeneous function of degree A (i.e. /(ax) = a^f{x.) for a > 0) and k :— 
— A + 1 — disa non-negative integer. This means / must decay with an exponent linear in the dimension 
d which does not apply to fields generated by a Riesz-kemel type source (see AppendixA.4). There are not 
only distributional contributions from the field derivatives, but from the fields themselves (some parts of 
the higher order field momenta). 



6.4 Ultraviolet-Analysis 

As for gauge theories in quantum field theory, dimensional regularization should be used in classical gen- 
eral relativity [31]. Therefore first all integrals must be evaluated in generic d dimensions and then the 
limit d ^ 3 is calculated. However, certain integrals are very difficult to solve for generic d. In practice 
one therefore evaluates the integrals in c? = 3 first and then determines possible additional contributions 
that arise from dimensional regularization. That is, one analyses the d-dependence of the integrals close 
to the singular sources, i.e., in the UV. (Only close to singularities regularization is important.) This is the 
purpose of the present section. Other necessary integration techniques are provided in Appendix A. 

The UV-analysis in generic dimension d is a necessary ingredient to correctly derive the Hamiltonians 
at formal 3PN level. This includes the 3PN point-mass Hamiltonian (see [31]) and the NNLO spin-orbit 
and spin(l)-spin(2) Hamiltonians considered in the present article. It would also be necessary for the yet 
unknown NNLO spin(l)-spin(l) Hamiltonian. 

For integrals only obtained for d = 3 one has no control on poles in — 3). There are two different 
problems with such poles: First the poles do not appear in pure d = 3 calculations and thus lead to 
ambiguous results after integrations by parts in integrands containing such poles (in one representation 
there are poles, in another maybe not). This comes from the fact that some of the pole terms can also give 
finite contributions which must be added to the d = 3 result. Second the poles have to cancel each other in 
order to extract a finite result from the d dimensional integration in the limit d — ?► 3 (or one must be able 
to absorb all poles through a renormalization procedure as in [184]). Both problems are well-known and 
also discussed in [31]. In the following we will provide some more technical details on how to perform the 
UV-analysis depending on the structure of the integrand. 

All integrals involving ft-^Jo) ij^ ^ki i^{2)^\^)) ™^ ^'S^ order potentials such as ^(g) or Tr^g^ are 
not available in d dimensions and were only calculated in d = 3 dimensions here. In all other integrals 
the limit d ^ 3 is straightforward, although integrations in d dimensions sometimes involve around one 
million terms on which the hmit must be performed. In case of the TT-projection of 0(2) '''(3) > ths fields 
are available in d dimensions. Hence, one can split up this part of the Hamiltonian in one-particle TT- 
projections (which can be performed in d dimensions) and two-particle TT-projections (which can only be 
evaluated for d = 3 with the presented methods). For the latter ones must still perform the UV-analysis. 

The term UV-analysis in this context refers to the short-range behavior of the integrand around a specific 
point. This will become more clear during the following explanation. Let us now consider the decay of the 
integrand f{ra,ri,, ria, rib) around the source a. First of all the integral is split up into a ball integral around 
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one of the sources, say for example source particle a, and an integral over the whole R'^ without this ball, 

/ d''x/(ra,r6,na,nf,) = / d'^a; /(ro, r;,, Ho, n;,) + / d"* a; /(ra, r;,, rio, n?,) , 

(126) 

where < £a ^ fab- The variables r?, and n;, of the other source (6 ^ a) are expressed in terms of r^, rab 
and Uab, 



rf, = |x - Zfcl = |x - Zq + Za - Zb| = y r2 + r;;, + 2rarab{'n.a ^ab) , (127) 

= —Ha H nab , (128) 

such that aU x-dependent expressions come from and ria type variables. Next we concentrate on the 
ball integral around a, 

/ dV(^a,n„) = /dr!,,rf_i / °dr,rf-V(r„n,). (129) 

JS£^(x„) J Jo 

Now the integrand is expanded in Va (leaving ria untouched). This is possible because a and b are well 
separated, and the ball contains only a small neighborhood of the source a. Then the integrand takes the 
form of a polynomial in ra and one can pick out the terms contributing poles at 3 (the ones with an 
exponent giving —3 for d = 3 on Tq). The next step is to count the number of riQ-vectors in each term and 
remove terms with an odd number of these vectors. This is due to the averaging procedure coming from the 
angular integration in (129) using the formulas (172) in the Appendix. Consider for example an integrand 
f{ra, Ha) = C {d)ra~^'^ {Pa nQ)(pfc ria) (this integrand is of the form qualified as dangerous in [31, Eq. 
(3.1)]) then (129) gives 
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= Cid)pa^PbJJ d^a,d-l« ^rarl~'"' . (130) 

Using (172b) and usual integration rules we obtain 

Cid) £2(3-d) 

A^xfira, no) = ^r!a,rf-i(Pa Pb) ^° . (131) 

The last integration step was performed by means of an analytic continuation from d < 3. At this stage 
there are three possibilities: The first one is that C{d) contains several factors of d — 3 which cancel the 
pole in the last factor and even lead to a vanishing limit when d — > 3. Then the potentially dangerous term 
is actually not dangerous at all. The second possibility is that C{d) ^ d ~ 3 which would also lead to a 
cancellation of the pole but would give a finite contribution which has to be taken into account to get the 
correct Hamiltonian. Last but not least C (d) could have such a structure that a pole remains and so this 
term has to be renormalized or canceled by another pole to give a physically meaningful result. 

The procedure mentioned above is only valid if there is no TT-projection (and in particular no ^(Jo)zj) 
appearing. The analysis of the hJJ^^^ ^ - type integrals works as follows (during this discussion we talk about 

a two-particle system where we consider ri as expansion point): The /i|7o) ij given in terms of 

inverse Laplacians and field variables in d dimensions by 

hf.o)., = -^^5ir^A-(0(,,,0(,),) , (132) 

see (122) for the explicit solution in d = 3. Now one can insert the 0(2) field, given by 
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where Ua = —IGnG'^'^^A ^Sa ^ After interchanging TT-projector and inverse Laplacian, one gets 

^IJ^'' i'^a,kUa,i) = 0. One can see this by using Ua.kUa.t ^ n^n^^i^^^^'^ which can be rewritten using 

d.d.ri-'^ = -2{d - 2)(<5,, - 2{d - l)nlni)rl-'^ , (134) 

as Ua.kUai ^ 2(d-i) i^ij^a^^'' + 2(d-2) ^'^J^a"^'^)' ^^^^ ^^^^ obviously bc projected to zero by the 
TT-projector Thus, after TT-projection there is only one part left, namely 

^(40) »j = -'^^^j-^mim2A-^6lJ'\ui^kU2A) ■ (135) 

Under the TT-projector one can integrate by parts because dkS^ — and obtains 



^(4o)*i = 2^|^^-4t"^i"^2A ^5lJ'^{uidkdiU2) . (136) 



Notice that the derivative should act on the U2 term because all quantities will be Taylor expanded in ri 
around x = Zi and ui is already proportional to an ri-power. At this stage of processing the hj^^-^ terms 
one is able to use the same analysis procedure as mentioned above: Check whether there are powers of ri 
with exponents smaller than —3 for d = 3 and expand dkdiU2 around ri to an order sufficient to reach the 
critical value in the exponent of ri . The idea is the expansion of the r2 variable in the TT-projector such that 
the Taylor expansion consists only of multiple inverse Laplacians on powers of ri, which can be calculated 
using (168) from the Appendix. Also bear in mind that the inverse Laplacian introduces additional powers 
of ri via the Green's function. The final ball integration can now be performed as discussed above. The 
same technique can also be used to perform the UV-analysis of terms involving S'^J'''' {(t>i^2)T^'^3))- 

To check our UV-analysis code we first reproduced the pole coefficients given in [31, Table 1]. The 
finite UV-contribution to the 3PN point-mass Hamiltonian in our representation (93) is given by 



.3PN.UV. ,^ _ 2A6-2'^(d - 2)(d + 1)(96 - 40d - 28rf^ + d^)Tr^~^^/^T (f ) 
™' 3(d-4)(d- l)4(d + 2) 



mim2 (d(ni2 pi)^ - P? + rf(ni2 P2)^ - Pa) , (137) 

''12 

where A is an UV-cutoff scale which does not contribute in the limit d 3. A similar analysis for the 
2PN point-mass Hamiltonian gave no contribution. We found no net contribution to the spin-dependent 
Hamiltonians, though poles and finite parts appeared in intermediate expressions. That is, Hadamard reg- 
ularization would have been sufficient to obtain the correct linear-in-spin Hamiltonians presented in the 
present work. The same situation was also found for the harmonic-gauge calculation of the equations of 
motion in [85, 86]. 

7 Results 

After discussing the several simplifications given above to reduce the integral of the formal 3PN Routhian 
to a form which can be handled appropriately, we continue by giving a short description of the integrands 
showing up at this order The integrands can be divided into three different types: 



the delta-type / d'^a;/(x)(5i, 
the Riesz-type J d'^x . . . n^-^nil . . . n^V^rj 
and the generalized Riesz-type j d'^x 



The solution of these three types of integrals will be shown in AppendixA. 

We used our Mathematica code to perform an integration of (93) directly with all fields inserted up 
to linear order in spin, neglecting spin(l)^ and spin(2)^ terms afterwards. From this we obtained the fully 
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reduced matter Hamiltonians for point masses, for the spin-orbit- and the spin(l)-spin(2)-interaction. There 
were no — 3) poles for the linear-in-spin Hamihonians giving rise to finite parts or being singular for 
d 3. Such poles only appeared in some intermediate steps in the UV-analysis (see Sect. 6.4) but finally 
identically canceled. For the point-mass parts there appeared a finite contribution given in (137). This 
together with our integration result reproduced the result from the literature exactly. 

Now we are able to write down the NNLO linear-in-spin fully reduced Hamiltonians in terms of matter 
variables only, by using (93) for the Routhian, inserting the sources (104), S(6)y (107b), Tr^attei (1 10)> ths 
solution for the constraint fields ( 1 1 1 )-( 1 1 5), (55), (56), and the propagating degrees of freedom (121) and 
(122). Notice that certain field derivatives may lead to distributional contributions mentioned in 6.3. The 
integrals can (for a binary) be solved by using the techniques given in AppendixA. Both Hamiltonians are 
valid for any compact objects like black holes or neutron stars. Their center-of-mass frame versions are 
given in 7.3 where the gyromagnetic ratios in the spin-orbit case are also given in [132, 133]. 



7.1 Next-to-next-to-leading Order Spin-Orbit Hamiltonian 

The spin-orbit Hamiltonian given in this subsection is the higher order gravitational analogue of the in- 
teraction of an electron's spin interacting with the electron's orbital angular momentum in the case of an 
e.g. hydrogen atom. In quantum electrodynamics this interaction is responsible for the fine structure in the 
spectrum. Here the spin is obviously no quantum mechanical quantity, it only characterizes the rotation 
(i.e. its direction and magnitude) of a gravitating mass in the gravitational field of another mass. Notice that 
the fine structure constant a in electromagnetic theory is substituted by Newton's gravitational constant G 
here. The result for the NNLO spin-orbit Hamiltonian reads 
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This Hamiltonian is formally at 3PN but for maximally rotating objects the post-Newtonian order goes up 
to 3.5PN. Recently in [85, 86] the NNLO spin-orbit contributions to the acceleration and spin-precession 
in harmonic gauge were calculated and agreement with the equations of motion following from our Hamil- 
tonian was found. ^ From a combinatorial point of view there are 66 algebraically different possible contri- 
butions to the Hamiltonian for each object (written in terms of the canonical spin tensor), but 24 of them 
do not appear in the canonical representation used here. 



7.2 Next-to-next-to-leading Order Spin(l)-Spin(2) Hamiltonian 

Also the spin(l)-spin(2) Hamiltonian has an electromagnetic counterpart. It is the gravitational analogue 
to e.g. the coupling between electron spin and spin of the atomic nucleus, responsible for the hyperfine 
structure in the electromagnetic spectrum. Of course in our case the spin(l)-spin(2) interaction leads to 
the modulation of gravitational waves but does not lead to a hyperfine structure in the emitted atomic 
electromagnetic spectrum. The result for the NNLO spin(l)-spin(2) Hamiltonian reads 
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This Hamiltonian is formally also at 3PN but for maximally rotating objects the post-Newtonian order goes 
up to 4PN. Notice that from a combinatorial point of view there are 167 algebraically different possible 
contributions to the Hamiltonian for all objects (written in terms of the canonical spin tensor), but 75 of 
them do not appear in the canonical representation used here. 

7.3 Hamiltonians in Center-Of-Mass Frame 

For later computations of, e.g., the mentioned orbital parametrizations of a binary system it is convenient 
to provide the Hamiltonians in the center-of-mass frame (pi = — P2 = p). In this frame in dimensionless 
quantities (see e.g. [96, 97] for rescaling) they are given by 
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There A = Si — S2 and S = Si + S2 the differences and sums of the spin vectors and L is the orbital 
angular momentum L = ri2ni2 x p. 



8 Kinematical Consistency: The Approximate Poincare Algebra 

For a space-time which is asymptotically flat the Poincare algebra must be fulfilled at spacial infinity, 
e.g. [189]. The generators of the Poincare algebra can be expressed in terms of the canonical variables 
describing the physical system, i.e. matter variables like linear momenta, position variables or spins. Also 
propagating field degrees of freedom enter the generators of the Poincare algebra. Throughout this section 
we set d — ?). The relations between the generators are given by 



{P„if} = 0, {J„il} = 0, 


(141a) 




(141b) 




(141c) 


{G„H} = P,, 


(141d) 


{G,,P,} = c-^d,,H, 


(141e) 


{Gi,Gj} = —c^'^eijkJk , 


(141f) 
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where P is the total Hnear momentum, is the total angular momentum tensor and Ji = ^eijkJ-'^ the 
associated dual vector, G is the center-of-mass vector and H the Hamiltonian of the physical system. Total 
linear momentum P and total angular momentum J*-' = — J^* are given by 



(142) 



see also, e.g., [30, 51]. For the contributions of the propagating field degrees of freedom see, e.g., [76, 77]. 
However, these contributions can be dropped within P and J*^ here as we are considering the conservative 
matter-only Hamiltonian instead of the ADM Hamiltonian (the latter still depends on the canonical field 
variables). 



8.1 General Considerations for a Center-Of-IMass Vector Ansatz 

As in [30, 51, 79] we use an ansatz for the center-of-mass vectors G at next-to-next-to-leading order (at 
lower orders it is also possible to directly calculate G from certain integrals). For constructing the center-of- 
mass vectors one has to consider the irreducible algebraic quantities which can be generated from Sa 
Pai and rj^f,. Since the Newtonian center-of-mass vector 

GN = ^maZa, (143) 



a a b^a 

is at the higher order corrections to the center-of-mass vector are also one post-Newtonian order below 
the appropriate Hamiltonian. Thus the momentum and G powers there are also reduced. 

Let us demonstrate these considerations at point-mass level: The Newtonian Hamiltonian has only 
terms at [which could be p^, pj, (ni2 Pi)^, (ni2 P2)^, and (pi P2)]^ and p° terms at G^ (which is 
only one term). At IPN there appear p"* terms at G°, terms at G^ and p'^ at G^. At 3PN there are p^ 
terms at and terms at G^. The center-of-mass vectors belonging to the Hamiltonians above have the 
following momentum powers emerging there: The Newtonian center-of-mass vector mentioned above has 
p° at G°, the IPN one contains at G° and p° at G^ and at 3PN level it contains p^ at G° and p° at G^. 

Now we will discuss how to construct the Unear-in-spin corrections for the center-of-mass vectors. 
Symbolically they can be written in the form 

Gso = SO-scalar • PM-vector + PM-scalar • SO-vector , (145) 
Gss = 5i52-scalar • PM-vector + SO-scalar • SO-vector + PM-scalai- • S'iS'2-vector . (146) 

Notice that we are formally working in generic dimension, where a spin-vector can not be defined. As the 
Poincare- Algebra must also hold in generic dimensions, it must be possible to construct the center-of-mass 
vector in terms of the spin-fensor. This is fortunate, as identities such as (6.1) in [85] would complicate the 
situation if one is forced to work with a spin vector in d = 3. 

Let us now summarize the vector quantities which can be built at certain spin-levels and may be used 
to construct the center-of-mass vectors. 

The mentioned vectors must perhaps be multiplied by scalar quantities, see (145) and (146). If the 
number of momentum variables in the spin-orbit or spin(l)-spin(2) scalars given in the following Table 2 
is not sufficient for the appropriate G-order they have to be filled up by the point-mass scalars, namely the 
linear momentum powers. 

Also important is that every spin is counted like a linear momentum because they have the same c^^- 
order, see 3.1. This means the formal 3PN spin-orbit and spin(l)-spin(2) Hamiltonians have only contri- 
butions up to G^ (G^ contributions cannot contain any spins since they are momentum-independent for 



^Although only pj and are appealing at the Newtonian order. This discussion should only provide some idea of appearing 
momentum powers at certain post-Newtonian orders. 
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Vector 


Irreducible Quantities 


point-mass (PM) 
spin-orbit (SO) (for Sa 
spin(o)-spin(&) (SaSb) 


Za, Pa, nab 

KbSa{i){j},PajSa(t)(j),PbjSait)U) 
Kb^a {k)U)Sb PakSa {k)U)Sb PbkSa (k)U)Sb 



Scalar 



PM 
SO 

S1S2 



Table 1 : Vector quantities at certain spin levels 
Irreducible Quantities 



linear momentum powers p" 

(KbPo. Sa (i) (j) ),^ {n\bPb j Sa (j) (j) )APa iPb j Sa (i) (j) ) 
{Pl tP2 J (j) Q) ) (pi iP2 j S2(i)U)), {n\2Pl J Si (,) (j-) )in\2PljS2(i)(J)), 
{n\2PljSj(i)(j)){n\2P2jS2(i)(j)), iPliP2jSl(i)(j)){piiP2jS2(i){j)), 

{nl^PijSi (i){j)){n\2PijS2 {n'i2Pi3Si {i){j)){n\2P2jS2 

{n\2P2 jSi (,) (j) ) {n\2Pi J S2 ), {n\2P2 jSi (,) q-) ) {n\2P2 j S2 (») (j) ), 

iSli^)(j)S2(^){j)), ("l2"i25'l(jK/c)'S'2(j)(fc)), 
{n\2Pl jSi (i)(fe)5'2(j)(fe)), (»^12P2j5'i(i)(fc)S'2(j)(fc)), 
{Plin{2Si (i)(fc)<5'2(j)(fe)), (P2ini2'^l {i)ik)S2 (i)(fc)), 
(PliP2i'S'l (i)(fe)52(j)(fe)), {p2zPljSi(i)(k)S2(j){k)),^ 
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in\2PljSi(^i)(^j)){piiP2jS2(i){j)), in\2P2jSi(i)(^j)){piiP2jS2(t)(j)) 

Table 2: Scalar quantities at certain spin levels 



point-masses). Notice that the Hamiltonians can only be constructed from the irreducible scalar quanti- 
ties given above. This is demanded by the Poincare algebra, namely (141a) (H should be invariant under 
translations and rotations and thus is a scalar). The Zq contribution in the center-of-mass vector can be 
fixed by (141e) using the lower order Hamiltonian. This Hamiltonian can always be written in the form 
H = J2a ^a, where the ha are translation invariant, {ha, P} — 0. (In the post-Newtonian approximation 
of general relativity all Hamiltonians have such a structure that the ha are translational invariant.) If we 
make an ansatz for the center-of-mass vector of the form 



G = ^ haia + Y , 

a 

we see that (c~^ — 1) 



(147) 



\ a b ) a 



{ha,P^}zl + haSabSij 
=0 



{Y\P^} = Y,haS,j=HS., 



=0 



(148) 



Equation (148) demands that {Y^,P^} = and so Y must be translational invariant. We have shown that 
the part of the center-of-mass vector which is not translation invariant, i.e., J^a ^a^a, can be read off from 
the Hamiltonian. 

From these consideration it follows that in the spin-orbit case the center-of-mass vector consists of 
52 algebraic independent quantities for one object and in the spin(l)-spin(2) case there are 86 algebraic 
independent quantities for both objects. Notice that up to the formal 3PN level and linear order in spin all 
center-of-mass vectors can be fixed uniquely by using the Poincare algebra. 



8.2 Next-to-next-to-leading Order linear-in-spin Center-Of-Mass Vectors 

Now we take the ansatz for the center-of-mass vector ( 147) where Y has to be constructed from the irre- 
ducible quantities given in Tables 1 and 2 with the decompositions (145) and (146), but without Za vectors 
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and put them into (141d) with the Hamiltonians (138) and (139) for the spin-orbit and spin(l)-spin(2) case. 
From this all unknown coefficients mentioned above could be fixed uniquely. The center-of-mass vector 
contributions given here implement the change in the binding energy of the system due to the NNLO linear- 
in-spin interaction Hamiltonians. This results also in context of the energy-mass equivalence in a modified 
gravitating mass and thus in a correction to the Newtonian center-of-mass vector Gn = J2a ''^a^a, which 
does not take any interactions into account. The correction to the center-of-mass vector from NNLO spin- 
orbit interactions finally results as 
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and the NNLO spin(l)-spin(2) part reads 
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From this the boost vector K. — G — tP can be obtained, which explicitly depends on time t. 

Notice that in (149) (in contrast to (150)) there appears a one-particle term without G factor It comes 
from the displacement of the center-of-mass due to the rotation and the resulting special relativistic Lorentz 
contractions of different parts of the object (which has to have a finite size). Further discussions about this 
issue can be found in [76] which are clarified graphically in particular in Fig. 1 therein. In (150) there is no 
term without G factor because all interactions between two spins are transmitted by the gravitational field 
in general relativity. 



9 Test-Spin near Kerr Black Hole 

In the last section we checked whether our results are compatible with kinematical restrictions of the 
Poincare algebra. Here we derive an exact test-spin Hamiltonian to compare our results (138) and (139) 
with. In the following subsections we restrict ourselves to rf = 3, because in the test-spin case there are only 
delta integrals to be evaluated. A partial check of our Hamiltonians against the test-spin case is contained 
in [74] for the case of aligned spins. 

There are various approaches to calculate the motion of a test-spin near a Kerr black hole (see, e.g., 
[190] and references therein). Since in the counting used in [ 190] the NNLO spin(l)-spin(2) interaction is 
at 4PN and therefore not considered therein, one needs to calculate the spin(l)-spin(2) contributions from 

ffxestspin - -i"e(,„)fee("] + J d^x [H'™""iV - nr''"N'] , (151) 

where N, Ni, the framefield e(„j)fe, and the implicit appearing metric provide the exterior gravitational 
field, and 7^™^'"=' and represent the test-spin moving in it [76]. Note that the framefield in the first 

term has to be evaluated at the position of the test-spin. 
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Since all spin dependencies of the metric are at least quadratic in the Kerr spin, and the only contribution 
which is linear in Kerr spin comes from the shift vector, the three-dimensional part of the metric and the 
lapse are identical with the appropriate components of the isotropic Schwarzschild metric which also comes 
from the spinless limit of the Kerr metric. The shift is given by the expressions in [60, Eq. (54)]. 

The metric components generated by particle '1' are given by 
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It is well-known that the metric components can be rewritten into a three-dimensional metric on the spatial 
hypersurface, lapse N, and shift — Nj, using (2) and (3). So one can immediately see that 
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with the squares of the shift neglected, because they are quadratic in the Kerr spin. 
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see [60]. Here ai (^^(j) = 5*1 /fni is the Kerr spin belonging to the black hole located at position '1'. 
Note that to linear order in spin 
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(159) 



(calculated from the inverse metric and the three-dimensional Christoffel symbols, see [60, Eq. (65)]) 
fulfills the ADM gauge condition, so no further coordinate shift from quasi isotropic coordinates to another 
coordinate system is necessary. Due to the symmetric framefield gauge e(i)j = .^TiJ, the framefield is 
given by (101) 
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For the test-spin in a Kerr field it is sufficient to calculate (151) where the sources are given by (94) and 
(95). Since the metric and the framefield are proportional to Sij many terms vanish in the source. The only 
terms remaining are 
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Inserting sources, metric components and framefield, and evaluating them at the testspin location gives 
the exact result 
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which leads after a post-Newtonian expansion (the post-Newtonian order given in the subscript is a formal 
one) 
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in full agreement with the test-spin limit of the full point-mass Hamiltonian, the spin-orbit Hamiltonian 
and the spin(l)-spin(2) Hamiltonian up to and including the formal 3PN order For later checks we provide 
also the expressions at formal 4PN order The test-spin limit for the next-to-next-to-next-to-leading order 
(NNNLO) spin-orbit interaction is given by 

Ken- Y45 mi(pi)3 , 13 mf(pi)^ , 315 mfpj ^ 105 mj V,,^ ^^^o^ 

4 

+ 14^((ni2 X p2)ai), (165) 

'^12 

and for the NNNLO spin(l)-spin(2) interaction by 
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''12 V / 

Notice that the formal 3PN spin(l)-spin(2) test-spin contributions were not given in [190] (in their counting 
rules they would be at 4PN level). Further notice that there are no contributions coming from the first term 
in (151) in the spin-orbit, and the spin(l)-spin(2) case, since in isotropic Schwarzschild coordinates it 
vanishes identically, see (96), (97), (155), (156), and (160). 

There are two possible further checks which use different approaches. As mentioned in [84] a com- 
parison of the effective field theory NNLO spin(l)-spin(2) potential in [82] to our NNLO spin(l)-spin(2) 
Hamiltonian would be a very strong check of both results since the EFT results are completely independent 
from the ADM formalism. Also a confirming check would be the derivation of both NNLO Hamiltonians 
using the spin-precession method shown in [51]. Due to the complicated structure of these comparisons 
they will be postponed to later publications. A very recent check of the NNLO spin-orbit Hamiltonian and 
the resulting equations of motion was performed in [85, 86] in harmonic gauge. Furthermore in [86] the 
near-zone metric was determined which is an important step towards the mentioned template calculations. 



10 Conclusions and Outlook 

We have derived the next-to-next-to-leading order spin-orbit and spin(l)-spin(2) Hamiltonians for binary 
systems. The spin-orbit Hamiltonian completes the knowledge of binary black hole dynamics up to and 
including 3.5PN order if the objects are rapidly rotating. For neutron stars also the leading order cubic-in- 
spin Hamiltonians are needed as the results in [60, 61] are valid for black holes only and tidal deformation 
effects become very important [71-73]. The Hamiltonians were checked using two methods. The fulfill- 
ment of the global approximate Poincare algebra was a major criterion for the correctness of the derived 
Hamiltonians in the extended ADM formalism. During this check the center-of-mass vectors could be de- 
termined uniquely from an ansatz. Since the approximate Poincare algebra is not sensitive to the static part 
of the spin(l)-spin(2) Hamiltonian and fixed only the difference of the two coefficients at the highest order 
in G of the spin-orbit Hamiltonian we performed further checks. The most simple test is a linear-in-spin 
approximation of the Hamiltonian of a test-spin moving near a stationary Kerr black hole. We rederived 
the test-spin Hamiltonian from [190] in a different manner in Sect. 9 (avoiding the use of Dirac brackets). 
A comparison was straightforward as the same gauge was used. A more elaborate test is the recalculation 
of both NNLO Hamiltonians via the spin-precession frequency method in [51] and will be part of a further 
publication. Also a comparison of the NNLO spin(l)-spin(2) Hamiltonian with the NNLO spin(l)-spin(2) 
potential given in [82] will be part of a further publication and would be a very strong check, because 
the derivation of this potential is completely independent from the ADM formalism. The most important 
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confirmation to date is the independent derivation of the NNLO spin-orbit acceleration in harmonic gauge 
[85]. 

The resuhs given in this article complete the knowledge of the post-Newtonian approximate dynamics 
for binary black holes up to 3.5PN. For general compact objects like neutron stars the leading order cubic- 
in-spin Hamiltonians are still unknown. The NNLO spin(l)-spin(2) Hamiltonian is at 4PN, if both objects 
are rapidly rotating, but there are still some tasks left to get the full post-Newtonian approximate dynamics 
up to and including 4PN. For general compact objects the leading order quartic-in-spin Hamiltonians are 
also unknown, they are only known for black holes. ^ Furthermore the NNLO spin(l)^ Hamiltonian - 
which is also at 4PN if the object is maximally rotating and maybe stronger than NNLO spin(l)-spin(2) 
- is completely unknown. Last but not least the up to the corrections to the 4PN point-mass 
Hamiltonian are also still unknown [191]. 

To get reasonable results for the templates the far-zone radiation field also must be calculated at higher 
order in the post-Newtonian approximation and also at higher orders in spin. The energy and angular mo- 
mentum loss are also not known at a post-Newtonian order corresponding to next-to-next-to-leading order 
linear-in-spin. If radiation and fluxes are known at such high orders also a parameterization is necessary. 
These three major ingredients are needed for the analytical description of gravitational wave templates 
which are very sensitive to higher order post-Newtonian and spin corrections. Analytical results are still 
important because for spinning binaries the parameter space (masses and spin-directions of the compo- 
nents) is very large and numerical simulations are so time consuming that they cannot be used to cover the 
whole parameter space. ^ 
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A Integration Techniques 



In this appendix we recapitulate the techniques needed to solve the integrals that appear during our calcula- 
tion. Most integration techniques are needed for d = 3 only, but whenever possible we provide results for 
generic d here as they may be used for checks. The short-range part of dimensional regularization using an 
UV-analysis is explained in Sect. 6.4. 



A.l Inverse Laplacians 

Inverse Laplacians are necessary to obtain solutions for the fields in the constraints, <j) and tt*. As one can 
see in (43), (45), and (36), the constraint equations always reduce to Poisson-like equations. On the right- 
hand-side of the equations there is always some source expression appearing. If the source only consists of 
Dirac deltas or their derivatives the inverse Laplacians can be calculated very easily, 

<>«) 

Also for some arbitrary power of Va one can calculate the inverse Laplacian immediately, namely 

j.a+2 

A"V" = . (168) 

The generalization to the application of multiple inverse Laplacians is straightforward. More interesting is 
the calculation of inverse Laplacians for powers of ra and which is a special case of the three-particle 
integral originally derived in [193, Eq. (7)] and also later given in [69, Eq. above (B14)]. The regular 
solution of A^^:p^ in d = 3 is well-known and given by 

A-^— ^\nSab^Hra + n + rab)- (169) 
ran 

It was first found in [194], also used in [32] and finally rederived in [69, Eq. (B14)]. (Notice that in 
[184, Appendix C] a d-dimensional generalization of (169) was given, but its extension to higher inverse 
Laplacians, i.e. calculation of A~"(r^^'^7'^^''), is highly non-trivial.) Such solutions are for example 
necessary for the derivation of hJJ^-^ [see (117)]. The inverse Laplacians for more general powers of ra 
and namely r™r^' {n,m> —1 and n, m odd) in d = 3 can be found by using the ansatz 

A-^rTr^ = Wr'''{ra,n,rab) + Wr'''{ra,rt,rab)lnSab, (170) 

where W^™'" and W^'"" are polynomials of (n+r7i+2)-th degree in ra, rb, and rab which consist altogether 
of 2( "'j'" +2)(m + n + 3) unknown coefficients. These coefficients have to be fixed by certain consistency 
conditions 

AA~^ = 1, AaA~^ - A"^Aa = 0, AfeA"^ - A-^Afe = 0, (171) 

where Aa — d[°'^ dl""^ denotes the Laplacian with respect to an object coordinate Za. These considerations 
can be generalized for higher inverse Laplacians. In [27] this technique was extensively used. (Generaliza- 
tions of the mentioned method are discussed in [195, Sect. V.A.] where the inverse Laplacians are denoted 
as superpotentials referring to a non-compact source. See also [35, Sects. 5.1., 5.2.] and references therein.) 

Before one of the inverse Laplacians discussed above can be applied, one has to get rid of possible 
ni and n2 vectors. A way to eliminate these n-vectors is to rewrite them as derivatives with respect to 
the particle coordinates, which are then commuted with the inverse Laplacian, see e.g. [196] and [78, 
Appendix C, Eqs. (C6)-(C8)]. 
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A.2 Averaging Procedures 

The averaging over the angle variables (necessary for Hadamard's finite part procedure and the local UV- 
analysis) is given by 



/ 



d^a,d-l = , (172a) 

y"dri,,d_i<<^ = ^na,d-lS^,i, , (172b) 

dfia,d-i<'<'<'<'' = ^(^^^_ 2) ^°-'^-^ (^»i»2'^»3»4 + Si.tjt^i^ + S,^,Ji^,^) , (172c) 



/ 



where (172d) for d = 3 coincides with Eq. (A28b) in [100]. r^a.d-i is the surface of the d-dimensional unit 
sphere around point a. We use the A^i-^^^ i^^ notation in the same manner like [100, Appendix A5]. This 
means ^{iii2-i«} ~ 12aee ^cr(ii)...(T(i£) where & is the smallest set of permutations {!,...,£) making 
^{iii2...ie} fully symmetrical in ii, . . . , ig. The factors in front of the deltas come from demanding that 
the trace over pairwise unit vectors should be one and the integral has the value of the surface of the unit 
sphere after tracing all unit vectors. An integral over an odd number of unit vectors is zero. 

Another necessary averaging procedure is the averaging in a subspace perpendicular to an axis given by 
a certain vector (this is necessary for reducing integrals of the Riesz-type and the generalized Riesz-type 
which may have a certain tensor structure to a linear combination of scalar integrals of the appropriate 
type). 

A d-dimensional vector space can be decomposed into a line characterized by a certain vector and a 
d — 1-dimensional subspace S being perpendicular to the vector Let this vector be a with = 1. Then 
the averaging of a tensor product of the vector m lying in the d — 1-dimensional subspace (its components 
parametrized by angular variables) should give only contributions if the number of vectors is even, namely 

(1>S-1, (173a) 
(m,)s=0, (173b) 
{■miTTij) s ^ apij , (173c) 



where {■■ ■)s denotes the d — 1-dimensional averaging and pij the projector onto the subspace. The 
projector fulfills the well-known identities 

Prj = Pjr , (174a) 
a%j^O, (174b) 
Pijpjk = Pik , (174c) 

and is given by pij — Sij ~ a^a^ . The identity in (173c) is motivated by realizing that {rairrij) s fulfills the 
first two projector identities. This can be rigorously shown under the conditions that m lies in the subspace 
perpendicular to a and |m| does not depend on the direction of m. If one contracts (173c) with (5y , then 
can be pulled out of the averaging brackets (as it is required to be constant in S) and thus 

Sij{mimj)s = = apu , (175) 

from which one can conclude that {pa — d — 1) 

a = . (176) 

d — 1 
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One can generalize the above result to 

. . . m^-)s = y^"^l;^?M (P^'^'^ ■ • • (177) 
(a — 3 + 2n)\\ 



A.3 Finite Part Integration 

In the Hamiltonian density there appear several integrals of the delta-type / d'^x /(x)(5a. The function /(x) 
is often singular at x = Za- In d = 3 these integrals are often evaluated by the well-known Hadamard's 
partie-finie (finite part) method (see e.g. [27, 32, 196]). Consider a function /(x) which is well-defined 
in a neighborhood of x = Zq and singular at this point. Then it is possible to expand this function in a 
Laurent series around this point using the auxiliary functions /„, 

oo 

/„(£) :=/(z„ + en) - ^ a„(n)e". (178) 

a=-N 

One defines the zero order coefficient of the Laurent expansion averaged over the n-vectors as the regular- 
ized value of / at Za, namely 

/reg(za) ■= J d^^aoin) . (179) 

See (172) for the appropriate averaging formulas. Notice that according to [29, 135] the finite part integra- 
tion is ambiguous at 3PN because in general 

(/l/2)reg(Za) ^ (/l )reg (Za ) (/2 )reg (Za ) , (180) 

consider e.g. ((/)(2)*)(za) vs. (0(2) (za))"* (see also [197, Appendix 2] for a discussion of this "tweedling of 
products" property). We will discuss how to avoid this issue in Sect. A. 5. 



A.4 The Riesz Kernel 

Because of the difficulties arising from using Dirac delta distributions as sources of fields in a non-linear 
theory like General Relativity (e.g. the product of distributions is not well-defined), reformulations of the 
delta distribution in terms of a function of finite width are also useful, e.g. the Riesz kernel. This kernel 
tends to a Dirac delta when the regulator (the finite width) tends to zero. The Riesz kernel is given by 

S''i-'-'^)- j2^l^f a-'- (181) 

Thus, in order to get rid of the singularities appearing due to the usage of Dirac deltas, one can replace 
them by the Riesz kernel and taking the limit — > for the Riesz kernel regulators after calculating the 
Hamiltonian. Another advantage of the Riesz kernel is that there will be no distributional contributions to 
some of the field derivatives mentioned in 6.3. Yet it is important to use Riesz kernels in generic dimension 
to avoid ambiguous results (as they appear in the finite part method). Unfortunately some of the integrals 
(in particular inverse Laplacians) for the solutions of the constraints or wave equation have not been solved 
to date if a Riesz kernel type source is used. Therefore we used the Riesz kernel method only when squares 
of delta distributions appeared or as a check for the other methods. Another problem of the Riesz kernel is 
that it breaks the general covariance of the theory explicitly via violating the contracted Bianchi identities 
y^jn^i/ _ g ^jjj Qjjg hopes that after the process — > the covariance is restored. This can be checked by 
using the Poincare algebra, see Sect. 8. See [198, above and below Eq. (94)] for a more detailed discussion 
of this issue. 
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A.5 Treating Contact Terms 



Integration of the delta-type integrands is a little subtle . It was mentioned in [3 1 ] that by analysis of Fourier 
representations of the integrals within dimensional regularization one can show that 

J d^xfif2hSa = (/l)reg(Za)(/2)reg(Za)(/3)reg(Za) , (182) 

up to and including the formal 3PN order (But to our knowledge no general proof was given yet.) We 
evaluated delta-type integrals with the help of this formula and used the d-dimensional Riesz kernel to 
calculate the regularized values of the functions at the source point appearing on the right hand side. We 
found that for all cases the regularized values of the fields agreed with results from the usual finite part 
method in d = 3, though this does not hold for products of the fields like {4>(^2)'^)i^a)- Further, the only 
field leading to nonvanishing finite parts in (93) is /i/Jqu, and its derivatives. 



A.6 Reduction to Scalar Integrals 

Besides the delta-type integrals discussed in the previous subsection (or after insertion of Riesz kernels), 
the only other type of integrals at formal 3PN level is of the form 

J d^x n\' ■ ■ ■ n\' ■ ■ ■ ■ n{' ■ r^r^sj^ . (183) 

The Si2 is introduced from derivatives of In si2, which in turn arise from certain inverse Laplacians. This 
subsection provides a formalism to reduce these integrands with complicated tensor structure to a linear 
combination of pure scalar integrands of the generalized Riesz-type (discussed in the next subsection). It 
is only possible to rewrite the vectors as derivatives with respect to the particle coordinates applied to a 
function of the type ^2 s]^2 (^^ for the inverse Laplacians) for the case 7 = 0,1,2..., because the S12 
also depends on the particle coordinate. For integer 7 = n > 

n k / \ / / \ 



Si2^{r,+r2+r,2T = l^l^\\{\ r{rr'r-^' , (184) 

fe=0 l=Q \ / \ / 

and therefore the mentioned functions can be reduced to a linear combination of products of the form 
''f ^2 ''12- This product structure is crucial to rewrite the vectors in (183) into particle derivatives. Thus 
one has to use another method to eliminate the vectors from the integrand for 7 < where the functions 
cannot reduced to a product form. In [27] a method is mentioned which can be used to get rid of angle in- 
tegrations resulting from the vectors by using an averaging procedure of the integrand in prolate spheroidal 
coordinates. We will present this procedure here in a slightly different way. 
First of all one has to get rid of, e.g., the n2 vectors by using the identity 

r\2 = — ni H ni2 . (185) 

7'2 f2 

Afterwards one still has to eliminate rii vectors from the integrand. By using rii — (x — zi)/ri and the 
orthogonal decomposition of rii with respect to ni2 one gets 



m = (m ni2) ni2 + , n^* = , (186) 



where 



{n^n^^)=:A = —^{rl^rl-rl^), (187) 



=: B 



2riri2 

— \l[{ri+r2Y-rl,] [rl, ^ {r^ - r^Y] . (188) 



2riri2 



and ~ Sij — n-i2"'i2 the projector onto the subspace perpendicular to ni2. Equation (187) can be 
derived by considering — — Z2 P — |x — zi + zi — Z2 P and expressing it in terms of ri, ri2 and A 
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itself. B is given by = n\'^'^^n\ = 1 — where we used the projector condition, n^^"^^ = 0, (186), 
and (187). Now a tensor product of such vectors can also be decomposed via (186) as 

d'^xr,,..,;,< = /d''xr,,...,,(n^+Ani2)*^...(n^+Ani2)'S (189) 



where i^. collects all constant tensors Uke linear momenta, spin tensors, and W\2 vectors. Notice that 
neither A nor B depend on the direction of \\.\ . Furthermore the integration extends over the whole 
such that one can rewrite the integrals of the type 

d''2:5,,..,,(ri,r2,ri2)n^^i...n^^'= , (190) 

in a more convenient way after a coordinate transform to prolate spheroidal coordinates. The Cartesian 
coordinates in terms of the prolate spheroidal coordinates 77, 0i , . . . , (j)d-2) are given by (generalization 
of the coordinates given in [199, pp. 752]) 

xi^fCv, (191a) 
- /v/(e'-l) (1-^7') cos ,^1, (191b) 



2^3 = fVi^^ - 1) (1 - V^) sin 01 cos 02 , (191c) 



Xd-l 



/v/(e'-l)(l-?7')sin0i...cos0rf_2 (191d) 



/V(e'-l)(l-?7')sin0i...sin0<i_2 (191e) 



where — l<77<l;l<(^<oo;0< 0d-2 < Stt; < <pi, . . . , 4>d~3 < . and 77 can be related to 
distance variables ri, ^2, ri2 by using 

. n +7-2 n -r2 ri2 

ri2 ri2 2 

Now one can easily perform the transformation of the volume element which is given by 

d'^x^ -7^^)d^d^dnd^2- (193) 

Thus we can perform an averaging over the d — 1-dimensional subspace perpendicular to ni2 by trans- 
forming the volume form 

d''xS.,....,{n,r2,n2)ni'^ ...ni'" = J dCd7^ fie - v')Sn...,iC,v) 
J dnd-2ni'' ...ni'^ , 

I dC d77 nd-2f{e - V^)Sn...^, V) (% ni ) , 

d'^x S,,...,,in,r2,n2) {ni'' . . .ni'") , (194) 



where {nj^"^) — and {n^'^n^-') = -^^'^^^ . Notice that has the same properties like p.^ and (. . . ) 
the same like (. . . )5 in A.2. Further notice that the average in the integrand only depends on ri, r2, and 
ri2. Thus one is able to reduce all integrals of the generalized Riesz-type containing n-vectors to a linear 
combination of scalar integrals. Notice that in the case 7 > the method from the present subsection is also 
advantageous over rewriting n-vectors as derivatives. The latter can lead to new singularities, logarithms, 
or even a need for a new kind of regulator, e.g., for an expression like 

d3xn»Sr„-3^ (195) 

By using the method of the present subsection these issues are avoided. 
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A.7 Integration using Generalized Riesz-Formula 

After reduction to scalar integrals one is left with integrals of the form J d'^x r"r2 s72- Po'" integer 7 > 
one can use (184) and end up with integrals whose solution were found by M. Riesz for arbitrary d 
[200, 201], namely 



2 l^a+P+d 

r(-f)r(-f)r(^i±i^) 



For d = 3 a generalization of the Riesz-formula was found by P. Jaranowski and G. Schafer in [27], 
. , /(o + 2TO.2)r(--a 4) ^ ^_ _^ _ ^ _ 

+ /i/2(/3 + 2,-/3-7-2) 

- /i/2(a + /3 + 4, -a - /3 - 7 - 4) - 1^ , (197) 

with 11/2(2;, y) := ^B(x^y)^ incomplete regularized Euler Beta function. One can express the in- 
complete Beta function (Euler integral of the first kind) i?i/2(a;, y) in terms of the GauB hypergeometric 
function 2F 1 

Bi/2{x,y) = ^2Fi + = ^i^i (x + y, x; 1 + -1) . (198) 

The function B{x,y) represents the Euler Beta function or Euler integral of the second kind with B{x, y) = 
^'-^^'^^ It turns out that the regularization procedure mentioned in [27, Eq. (B21), (B23), and (B24)] only 



r{x)r(y)- 

modifies a and /3 to be non-integer via the introduced analytical regulators ^le and i^e. So one can in 
principle simplify the formula given above by the assumption that 7 e Z, and a and f3 being arbitrary. 
Positive integer powers of S12 can be handled through (184) and (196). Thus the only relevant powers of 
S12 are the negative integer ones. 

Equation (198) leads to hypergeometric functions of the type 2fi(— 7, z;z + l; —1) in (197), where z 
depends on a and /3, namely 

d^xr'^r^sj, = 2nr{2 + a)T{2 + /?)r(-4 - a - (3 - 7)|-^^ 

2j^i(-7,2 + a;3 + a;-l) 2^1 (-7, 2 + /?; 3 + /3; -l) 
r(3-Ka)r(-2-a-7) r(3 + /3)r(-2 - /3 - 7) 

_ 2^1 (-7, 4 + g + /?; 5 + g + ~1) I ^^+^+^+3 .j^g. 
r(5 + a + /3)r(-4-a-/3-7) / 12 • 

So for negative integer gamma 7 = — n one can express the solution of the integral in terms of 2F1 {n,z; z + 
1; -1). It is well-known that [202] 

2Fi(0,z;z + l;-l) = 1, (200) 

z 



2Fi(l,z;z + l;-l) = ^ 



(201) 



where ip is the Digamma function. From these two formulas and the contiguous relation of the GauB 
hypergeometric function [199, 203, 204] 

= (c — 0)2-^1(0 — 1, b; c; z) + (2a — c + (6 — 0)2)2-^1(0, b; c; z) 
+ o(z- 1)2^1(0+ 1,6; c;z), (202) 
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one obtains the recursive relation 



2Fi(n,z;z + l;-l) = 



(z - (n - 2))2Fi(n - 2, z; z + I; -I) 



2{n- 1) 

+ (3n - 2(2 + z))2Fi (n - 1, z; z + 1; -1) 



(203) 



From these relations one can show per induction the general structure of the hypergeometric function to be 



2Fi{n,z;z + 1;-1) = a}! 



(-ir 



k=Q 



Unfortunately the coefficients a^"^ can only be obtained by complicated recursive relations between differ- 
ent k and n following from (203) and we can only give an explicit form for some of them 



2(n- 1)! 



z + 1 



-1 IL— 1 

^(1) n(^-^)- (204) 



k=0 



An) _ 



An) 



0,n>2 and4°^ = 1 ,4^' = 0, 
(-1)" 



, n > 2 . 



2(n-l)! 
(„) _ (-l)"(l + n-n^) 
'""2" 4(n-l)! 



n > 3. 



(205) 
(206) 

(207) 



Nevertheless the mentioned recursion relations can be used to cache all GauB hypergeometric functions. 
Mathematica is able to handle limits, series expansions and derivatives of the arising Digamma functions 
very well and thus all occurring integrals in the binary Hamiltonian to linear order in the spin variables can 
be solved. 
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